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1.0  ABSTRACT 


This  report  describes  an  investiqation  into  the  problem  of  the  "exact" 
calculation  of  three-dimensional  lifting  potential  flows.  The  designation 
"exact"  is  used  to  denote  a  method  that  makes  no  approximations  in  its  basic 
formulation,  such  as  small-perturbation  or  lifting-surface  theories  do. 
Obviously,  numerical  realities  require  some  approximate  techniques  in  the 
computer,  but  "exact"  methods  can  be  numerically  refined  in  principle  to  give 
any  degree  of  accuracy. 

The  first  part  of  the  study  is  a  look  at  the  problem  of  three-dimensional 
lifting  potential  flow  from  a  fundamental  standpoint,  something  almost  totally 
lacking  in  the  literature.  Unlike  nonlifting  flow  whose  "physics"  and  mathe¬ 
matical  description  seem  basically  related,  the  mathematical  description  of 
the  lifting  problem  is  merely  a  model  to  describe  by  means  of  an  inviscid  flow 
a  phenomenon  that  is  ultimately  due  to  viscosity.  This  Is  true  even  in  two 
dimensions,  but  in  three  dimensions  It  leads  to  certain  logical  difficulties. 

The  method  of  this  report  and  all  current  "exact"  methods  of  calculating 
lifting  flows  are  based  on  the  author's  previous  work  on  three-dimensional 
nonlifting  flows.  This  report  describes  the  present  method  in  general  and 
In  detail,  including  all  formulas  and  logic.  Alternatives  are  discussed,  some 
of  which  are  discarded,  while  others  are  Incorporated  Into  the  program.  The 
present  method  differs  from  other  current  methods  mainly  in  its  use  of  finite- 
strength  surface  vorticity  distributions  instead  of  concentrated  line  vorticity 
interior  to  the  body  and  in  its  application  of  the  Kutta  condition.  Comparisons 
indicate  advantages  for  the  formulation  of  the  present  method. 

A  variety  of  cases  calculated  by  the  present  method  are  presented  to 
Illustrate  its  versatility  and  usefulness.  Comparisons  of  the  calculations 
with  experimental  data  are  presented.  The  Importance  of  viscosity  in  the 
experimental  results  is  illustrated. 
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4.0  PRINCIPAL  NOTATION 
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ij 
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C 
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P 


c 


c 


z 


d 

F,  S 


I 


nm 


i 


ij 


velocity  induced  at  the  i-th  control  point  by  a  unit  value  of 
source  density  on  the  j-th  element.  If  there  are  N  on-body 
control  points  where  a  normal -velocity  boundary  condition  is 
applied,  this  is  an  N  x  N  matrix.  It  is  the  coefficient 
matrix  for  the  linear  equations  for  the  values  of  source 
density.  The  same  coefficient  matrix  applies  to  all  onset 
flows.  j 

Constant  of  proportionality  for  the  dipole  strength  along  an 
N-line.  Local  dipole  strength  along  an  N-line  equals  B  times 
the  arc  length  along  the  N-line  from  the  trailing  edge.  By 
theorem  of  Appendix  A,  this  means  B  equals  the  value  of 
bound  vorticity  at  the  spanwise  location  of  the  N-line.  Used 
with  superscript  k  to  Indicate  value  {of  3  at  the  midspan 
of  the  k-th  lifting  strip.  j 

intercepts  of  slanted  sides  of  a  trapezoidal  element  with  the 
x-axis  of  its  own  coordinate  system  (figure  20). 

lift  coefficient  for  a  complete  body.  ^ 

pressure  coefficient.  Equals  difference  of  local  static  pres¬ 
sure  from  freestream  static  pressure  divided  by  freestream 
dynamic  pressure. 

denotes  an  integration  path.  Also  a  constant  multiplying  a 
second  order  dipole  term  used  to  produce  continuity. 

section  lift  coefficient.  Lift  force  on  a  strip  of  elements 
on  a  wing  divided  by  the  projected  area  or  thb  strip  in  a 
olane  containing  the  chord  line  and  by  freestream  lynamic 
pressure. 

used  with  double  subscript  to  denote  length  of  a  side  of  a 
quadrilateral  element.  ; 

subscripts  and  superscripts  used  to  denote  quantities  associ¬ 
ated  with  the  two  N-lines  bounding  a  strip  of  elements.  F 
denotes  "first"  N-line  and  5  the  "second"iN-\i ne .  f  is 
also  used  to  denote  number  of  uniform  onset! ■ lows . 

normalized  moment  of  the  arrc  nf  a  tra.r-v^AU'.  elem  wit! 
respect  to  the  axis  of  the  coordinate:  equa¬ 

tions  (7.2.24)  and 

a  subscript  used  to  denote  guaot.i *  ies  associate*'.  *.i i:n  the 

i-th  control  point,  partioul  a'ly  ve’ac Sties  at  the  t.  point. 

Used  as  superscript  to  denote  in, put  point. 

double  subscript  used  to  rienota  of  }«+.*»  o'- • 

i-th  control  point,  particular  *y  irvh.r./r.  velt.r-’ty. 
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i~£,  Tr»  Tc  unit  vectors  along  the  axes  of  a  coordinate  system  based  on 
an  element. 

k  subscript  used  in  various  ways,  k  =  1,  2 , 1 3 ,  4  denotes 

quantities  associated  with  the  fouij  corner  points  of  an  element. 
Also  used  as  subscript  and  superscript  to  denote  k-th  lifting 
strip  or  vorticity  onset  flow  associated  with  that  strip. 

L  arc  length  along  an  N-line.  Also  denotes  iota!  number  of 

lifting  strips  of  surface  elements.-  ! 

M  used  in  figure  42  to  denote  freestream  Mach  number 

m32’  m41  slopes  of  the  slanted  sides  of  a  trapezoidal  element  with 

^  '  respect  to  the  y-axis  of  its  own  coordinate  system  (figure  20). 

N  total  number  of  surface  elements  at  which  normal -velocity  boundary 

conditions  are  applied.  Includes  both  lifting  and  nonlifting 
elements. 

N-llne  curve  in  wing  surface,  usually  a  fixed  spanwise  location,  along 

which  input  points  are  given.  N-line  continues  aft  to  define 
the  trailing  vortex  wake.  A  strip  of  elements  lies  between 
two  consecutive  N-lines. 

n"  unit  normal  vector.  * 

0  number  of  off-body  points  at  which  flow  is  to  be  computed 

P  a  general  point  In  space. 

r  distance  between  two  points.  Used  with  subscript  o  to  denote 

distance  from  centroid  of  an  element  to  point  where  velocity 
Is  being  computed.  Used  with  subscript  k  to  denote  distance 
between  such  a  point  and  the  corner  point  of  an  element. 

S  denotes  a  body  surface  on  which  a  normal -velocity  boundary 

condition  is  applied. 

s  arc  length,  especially  arc  length  along  ar  N-line. 

:num  diagonal  of  an  element  (figure  20)1. 
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perturbation  velocity  due  to  body. 

total  flow  velocity  at  i-th  control  point  due  to  flow  induced 
about  the  body  by  the  k-th  vorticity  onset  flov..  With  super¬ 
script  (®)  the  nonlifting  flow  about  the  body  in  a  uniform 
freestream,equation  (7.13.4). 

velocity  components  induced  by  an  element  at  a  point  space 
with  respect  to  the  coordinate  system  of  the  element. 

velocity  Induced  at  the  i-th  control  point  by  a  unit  value  of 
source  density  on  the  j-th  element. 

velocities  induced  at  the  i-th  control  point  by  a  dipole 
distribution  on  the  j-th  trapezoidal  element  that  varies 
linearally  from  zero  on  one  parallel  side  to  unity  on  the 
other.  Superscript  denotes  the  N-line  containing  the  side 
with  nonzero  dipole  strength. 

width  of  a  trapezoidal  element  in  direction  normal  to  the 
parallel  sides  (figure  20).  Also  used  with  subscript  k  to 
denote  width  of  lifting  strip  for  parabolic  fit  (section  7.11). 

coordinates  of  a  point  in  element  coordinate  system. 

coordinates  of  a  point  in  the  reference  coordinate  system  used 
to  Input  the  body. 

coordinates  of  the  centroid  of  an  element  in  the  reference 
coordinate  system. 

direction  cosines  of  a  point  in  space  with  respect  to  the 
coordinate  system  of  an  element  based  on  the  centroid  as 
origin.  Also  used  with  subscript  k  to  denote  the  same 
direction  cosines  with  origin  shifted  to  a  corner  point. 

total  circulation  around  a  closed  path. 

circulation  about  a  c.osed  path  due  to  perturbation  velocity 
field  of  the  body. 

dipole  strength  per  unit  area. 

x,  y  coordinates  of  a  point  of  an  element  in  Its  own  coordinate 
system.  Used  with  subscripts  k  to  denote  coordinates  of  the 
corner  r-elnts. 

o stance  criteria  used  to  decide  when  multipole  and  far-field 
formulas  are  to  be  used. 


o  source  density  per  unit  area.  Used  with  subscript  j  to  denote 

value  on  j-th  element  and  with  superscript  k  to  denote  values 
calculated  for  k-th  vorticity  onset  flow. 
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velocity  potential  especially  that  due  to  a  body  or  that  due 
to  a  surface  element. 

velocity  potential  due  to  a  dipole  distribution  on  an  element 
that  varies  as  the  p-th  power  of  5  and  the  q-th  power  of 
rij  equation  (7,4.4). 
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5.0  INTRODUCTION 


5.1  Statement  of  the  Problem  of  Potential  Flow 

The  problem  considered  Is  that  of  the  flow  of  an  Incompressible  inviscid 
fluid  in  the  region  R'  exterior  to  (or  interior  to)  a  given  boundary  surface 
S.  For  definiteness  S  is  shown  as  a  single  three-dimensional  surface  in 
figure  1,  but  S  may  consist  of  several  disjoint  surfaces,  and  the  problem 
may  be  either  two-  or  three-dimensional.  It  is  convenient  to  express  the 
fluid  velocity  field  V*  at  any  point  P  as  the  sum  of  two  velocities: 

?  =  V  +  7  (5.1.1) 

oo 

The  velocity  is  denoted  the  onset  flow  and  is  defined  as  the  velocity 
field  that  would  exist  if  all  boundaries  were  simply  transparent  to  fluid 
motion.  It  is  assumed  that  V  is  known.  Most  commonly  V  represents  a 

oo  oo 

uniform  parallel  stream  and  is  thus  a  constant  vector.  The  vector  \T  is  the 
disturbance  velocity  field  due  to  the  boundary  surface  S.  Since  the  flow  is 
incompressible,  both  V  and  v  have  zero  divergence.  It  is  further  assumed 


\ 


Figure  1.  Potential  flow  about  a  three-dimensional  body. 
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that  v  is  Irrotatlonal ,  l.e.,  has  zero  curl.  Thus,  v  may  be  expressed  as 
the  negative  gradient  of  a  potential  function  $>, 

v  =  -grad  *  (5.1.2) 

The  con  Jon  of  zero  divergence  then  yields  Laplace's  equation  for  <j>, 

v24,  =  0  in  R*  (5.1.3) 

The  boundary  condition  on  S  is  derived  from  the  requirement  that  on  a 
stationary  Impervious  surface  S  the  normal  component  of  fluid  velocity  must 
vanish.  Thus, 


1^-  =  grad  $  •  rT  =  V"  •  r?  on  S  (5.1.4) 

an  “ 

where  rT  Is  the  unit  outward  normal  vtctor  to  S.  Since  the  right  side  is 
known,  equation  (5.1.4)  expresses  a  Neumann  boundary  condition  for  If  the 
boundary  S  is  moving  or  if  a  nonzero  normal  velocity  Is  prescribed,  the 
right  side  of  (5.1.4)  is  modified  in  an  obvious  way. 


A  regularity  condition  at  infinity  is  also  required.  In  the  usual 
exterior  problem  the  condition  is 

| grad  4>|  0  at  Infinity  (5.1.5) 

In  addition  to  the  above  equations,  some  applications  require  certain  auxil¬ 
iary  conditions  to  be  satisfied.  However,  In  the  absence  of  such  conditions 
and  for  a  simply  connected  region  R',  the  equations  (5.1.3),  (5.1.4),  and 
(5.1.5)  comprise  a  well-posed  problem  for  the  potential  <p. 

In  two-dimensional  exterior  problems,  the  region  R'  is  not  simply 
connected,  and  equations  (5.1.2),  (5.1.3),  (5.1.4),  and  (5.1.5)  do  not  define 
a  unique  velocity  field.  Define  the  total  circulation  r  around  any  closed 
path  c  in  the  fluid  as  the  line  integral 

r  =  JV  •  Hs *  =  •  3s"  +  J'v  •  ds  rm  +  y  (5.1.6) 

c  c  c 
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where 


(5.1.7) 


is  the  circulation  associated  with  the  disturbance  velocity  due  to  the  body. 
In  the  above 


ds  =  t  ds  (5.1.8) 

where  S  Is  arc  length  along  c,  and  t  Is  the  unit  tangent  vector.  If  c 
does  not  enclose  all  or  part  of  S,  then  y  =0.  If  S  is  a  single  surface, 
it  can  be  shown  (reference  1)  that  the  velocity  field  v  is  rendered  unique 
by  specifying  y  for  any  c  that  encloses  S.  If  S  consists  of  several 
disjoint  surfaces,  y  must  be  specified  for  a  set  of  paths,  each  of  which 
encloses  exactly  one  of  the  disjoint  surfaces  that  comprise  S.  The  potential 
4>  is  unique  If  and  only  if  y  =  0  for  all  closed  paths. 

5.2  Potential-Flow  Model  for  Lift 

The  reasoning  leading  up  to  the  formulation  of  the  potential  flow  problem 
in  terms  of  equations  (5.1.3),  (5.1.4),  and  (5.1.5)  seems  very  plausible. 

However,  when  the  problem  defined  by  these  equations  is  solved,  the  resulting 
flow  gives  zero  net  force  on  a  closed  three-dimensional  body.  This  is  due  to 
the  fact  that  all  components  of  force  cn  a  body  -  both  the  lift,  which  is  per¬ 
pendicular  to  the  freestream,  and  the  drag,  which  is  parallel  to  the  freestream  - 
are  ultimately  due  to  viscosity.  Nevertheless,  the  goal  of  calculating  at  least 
the  lift  component  of  the  force  by  a  purely  inviscid  technique  has  been  con¬ 
tinuously  pursued.  It  Is  important  to  realize  that  any  such  formulation  is 
simply  a  potential-flow  model  of  real  lifting  flow,  and  that  the  two  flows 
are  not  necessarily  related  in  any  fundamental  way.  Formulation  of  the  commonly 
accepted  potential-flow  model  of  three-dimensional  lifting  flow  has  relied 
heavily  on  results  for  the  two-dimensional  case. 

In  two-dimensional  flow  advantage  can  be  taken  of  the  indeterminacy 
of  the  solution  as  described  in  section  5.1.  For  a  single  closed  body  in  a 
uniform  stream,  the  drag  force  is  zero,  and  the  lift  is  proportional  to  the 
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circulation  y,  whicn  is  arbitrary.  (For  a  uniform  onset  flow  the  total  cir¬ 
culation  r  equals  y,  the  circulation  due  to  the  disturbance  velocity.) 

Thus,  in  two-dimensions  the  problem  is  not  that  no  lift  is  obtained  but  that 
the  lift  can  have  ai.y  magnitude.  Some  auxiliary  condition  is  necked  to  fix 
the  value  of  lift.  For  bodies  with  continuous  slope  no  satisfactory  auxiliary 
condition  has  ever  been  formulated.  However,  a  conventional  airfoil  has  a 
sharp  corner  at  its  trailing  edge,  and  there  is  a  unique  value  of  y  (and  thus 
a  unique  lift)  that  makes  the  potential-flow  surface  velocity  finite  at  this 
corner.  Determining  the  value  of  circulation  in  this  way  also  insures  that  a 
streamline  of  the  flow  leaves  the  airfoil  at  the  trailing  edge  with  a  direction 
along  the  bisecto*'  of  the  trailing-edge.  This  condition  of  finite  velocity 
at  the  trailing  edge,  the  so-called  Kutta  condition,  is  so  well  accepted  that 
it  is  normally  not  considered  a  mere  modeling  device  but  is  assumed  to  have 
a  more  fundamental  connection  with  the  real  flow.  However,  the  Kutta  condi¬ 
tion  is  inapplicable  to  smooth  bodies,  and  for  airfoils  with  sharp  trailing 
edges  it  gives  values  of  lift  that  differ  from  experimental  values  by  up  to 
20  percent. 

The  theorem  that  guarantees  a  unique  solution  for  the  flow  about  a  two- 
dimensional  body  with  prescribed  circulation  y  is  quite  general.  However, 
in  a  specific  calculation  procedure  the  question  arises  of  how  the  condition 
of  prescribed  circulation  Is  to  be  applied.  All  procedures  accomplish  this 
with  the  help  of  vortlcity.  A  distribution  of  vortlclty,  consisting  of  either 
concentrated  filaments  or  finite-strength  surface  or  volume  distributions  are 
hypothesized  to  lie  on  or  within  the  body  In  question.  The  total  strength  of 
the  vortlcity  distribution  establishes  the  prescribed  circulation. 

Consideration  of  the  above  two-dimensional  model  suggests  certain  elements 
of  a  model  for  lifting  flow  about  a  three-dimensional  wing  of  the  type  shown 
in  figure  2.  If  the  trailing  edge  of  the  wing  is  a  sharp  corner,  a  plausible 
three-dimensional  Kutta  condition  requires  that  the  velocity  remain  finite 
there  all  across  the  span,  which  means  that  a  stream  surface  leaves  the 
wing  from  the  trailing  edge.  Define  the  circulation  about  a  particular  wing 
section  as  the  line  integral  of  the  velocity  in  the  form  of  equation  (5.1.7) 
about  a  closed  curve  lying  In  the  wing  surface  as  shown  in  figure  2.  The 
precise  definition  of  this  so-called  section  curve  is  not  considered  now.  A 
reasonable  definition  is  that  the  curve  lie  in  a  plane  parallel  to  the  plane 


SPANWlSE 

direction 
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Figure  2.  Nomenclature  for  a  three-dimensional  wing. 

of  symmetry  of  the  wing.  But  for  certain  purposes  the  curve  could  lie  in  a 
plane  normal  to  the  leading  or  trailing  edge.  In  any  case  the  value  of  the 
circulation  is  different  for  curves  at  different  locations,  so  that  there  is 
a  "spanwlse"  variation  In  "section  circulation."  By  analogy  with  two- 
dimensions,  It  Is  expected  that  a  proper  adjustment  of  this  spanwise  variation 
could  render  the  velocity  finite  all  along  the  trailing  edge.  Presumably,  the 
circulation  can  be  generated  by  some  distribution  of  vorticity  lying  on  or 
within  the  wing.  It  seems  evident  that  the  direction  of  this  so-called 
"bound  vorticity"  should  be  generally  along  the  span,  roughly  parallel  to  the 
trailing  edge.  The  net  vorticity  strength  through  each  "section"  is 
proportional  to  the  circulation  around  that  section. 

Define  y-j  and  ^  as  the  values  of  circulation  about  two  sections  of 
the  wing,  where  the  positive  sense  of  the  integral  of  (5.1,7)  is  taken  as 
clockwise  to  an  observer  at  the  wing  midplane  looking  towards  the  right  wing 
tip.  Unlike  the  two-dimensional  case,  the  region  exterior  to  a  closed  three- 
dimensional  body  is  simply  connected,  so  that  if  the  flow  is  potential,  i.e., 
has  zero  curl,  and  is  free  from  singularities,  then 


c 
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for  any  closed  path  c,  which  implies  yj  =  yg  =  0.  Thus,  to  obtain  nonzero 
values  of  section  circulation,  there  must  be  some  form  of  singularity  in  the 
exterior  flow.  The  nature  of  the  singularity  can  be  exhibited  by  considering 
the  path  c  shown  in  figure  3a.  The  line  integral  of  velocity  around  this 
path  is 


J  V  •  <Js  =  Yj  -  Yj,  +  J"  ( v+  -  v_ )  •  37  (5.2.2) 

c  I 

where  I  is  the  straight  path  joining  the  two  section  curves  and  7+  and  7 
are  the  limiting  velocities  obtained  by  approaching  I  from  two  different 
directions  on  the  surface.  If  the  line  integral  of  (5.2.2)  is  to  vanish,  then 
either  yj  =  Y^  or  v*+  t  v"  ,  and  there  Is  a  discontinuity  of  tangential 
velocity  along  I.  If  sharp  corners  in  streamlines  are  to  be  avoided,  such  a 
discontinuity  can  occur  only  across  a  stream  surface  of  the  flow,  and  thus 
either  I  is  a  locus  from  which  a  stream  surface  leaves  or  joins  the  body  or 
else  I  Is  a  portion  of  a  streamline  on  the  surface.  In  any  event  I  repre¬ 
sents  the  intersection  of  a  sheet  of  vorticity  with  the  body  surface.  To 
complete  the  potential  flow  model,  the  first  possibility,  a  stream  surface 


(a)  (b) 


Figure  3.  Circulation  on  a  three-dimensional  wing,  (a)  Integration 
path  c.  (b)  Discontinuity  at  the  trailing  edge. 
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leaving  the  body,  Is  selected,  essentially  on  physical  grounds.  It  Is 
reasoned  that  vortlclty  Is  Introduced  only  to  the  fluid  that  passes  by  the  body 
and  that  the  path  I  of  (5.2.2)  must  lie  along  the  trailing  edge  of  the  wing 
(figure  3b).  Thus,  a  vortex  sheet  issues  from  the  trailing  edge  and  for 
steady  flow  it  proceeds  to  infinity.  The  average  strength  of  the  sheet  along 
I  is  proportional  to  the  difference  y-|  “  >2*  taking  the  limit  as  the  two 
section  curves  approach  each  other  gives  the  result  that  the  local  strength 
of  the  trailing  vortex  sheet  is  proportional  to  the  "spanwise"  derivative  of 
the  "section  circulation." 

It  follows  from  the  above  that  the  local  strength  of  the  "trailing 
vortlcity"  that  issues  from  the  wing  trailing  edge  equals  the  "spanwise" 
derivative  of  the  "bound  vorticity."  Thus,  trailing  vortlcity  is  of  precisely 
the  right  form  so  that  the  entire  oound-plus-trailing  vorticity  system  may 
be  thought  of  as  being  composed  of  constant-strength  vortex  lines  of  infin¬ 
itesimal  strength,  each  of  which  proceeds  "spanwise"  along  the  wing  and  then 
turns  and  proceeds  "streamwise"  to  infinity,  the  familiar  "horseshoe"  vortices. 
This  is  crucial  because,  as  pointed  out  in  reference  2,  the  velocity  field 
due  to  a  variable-strength  vortex  filament  or  a  nonclosed  constant-strength 
vortex  filament  of  finite  length  is  not  a  potential  flow.  Only  infinite  or 
closed  vortex  lines  of  constant  strength  give  rise  to  irrotational  velocity 
fields. 

As  mentioned  above,  the  trailing  vortex  sheet  must  he  a  stream  surface  of 
the  flow.  Also,  on  physical  grounds  the  pressure  must  be  continuous  across 
the  sheet.  In  principle,  these  two  conditions  allow  the  complete  shape  of  the 
trailing  vortex  sheet  to  be  calculated.  The  basic  flow  problem  is  nonlinear 
because  the  location  of  the  sheet  changes  for  different  onset  flows.  In 
particular,  the  sheet  changes  location  if  the  angle  of  attack  of  the  freestream 
changes. 

The  above  contains  the  general  features  of  the  potential-flow  model  of 
three-dimensional  lift.  It  is  considerably  more  complicated  than  the  simple 
formulation  of  eguations  (5.1.3),  (5.1.4),  and  (5.1.5),  which  represent  the 
nonlifting  case.  However,  the  nonlifting  formulation  appears  to  be  fundamental, 
while  the  lifting  formulation  is  basically  a  model  adopted  to  simulate  certain 
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properties  of  real  viscous  flow  by  means  of  a  potential  flow.  The  nonfunda- 
mental  nature  of  the  lifting  model  leads  to  some  logical  difficulties  which 
may  or  may  not  be  Important  In  a  particular  case.  Some  of  these  are  discussed 
In  the  next  section. 

5.3  Some  Logical  Difficulties  in  the  Potential-Flow  Model 

The  principal  device  by  which  lift  is  introduced  into  potential  flow  of 
either  two  or  three  dimensions  is  the  trailing  edge.  To  some  extent  the  defini¬ 
tion  of  a  trailing  edge  Is  a  matter  of  legislation  by  the  user  of  the  method 
rather  than  a  fundamental  concept.  Accordingly,  difficulties  may  arise.  In 
two-dimensions  the  situation  is  rather  simple.  There  is  no  logical  difficulty 
if  the  trailing  edge  is  a  sharp  corner  (the  agreement  of  the  model  with  real 
flow  may  or  may  not  be  acceptable).  On  the  other  hand,  if  there  is  no  sharp  ] 

corner,  the  difficulty  is  crucial,  because  the  trailing  edge  cannot  be  1 

rationally  defined.  In  three-dimensions  some  rather  subtle  borderline  cases  | 

arise  in  ordinary  design  applications.  In  regions  where  the  wing  has  a  sharp  j 

corner  as  shown  in  figure  2,  the  choice  of  trailing  edge  is  straightforward.  | 

Difficulty  arises  where  the  locus  of  the  sharp  corner  ends.  The  question  ' 

arises  whether  the  trailing  edge  ends  or  continues,  and,  if  the  latter,  in 
what  matter.  I 

' 

A  wing  tip  Is  the  place  where  the  above-mentioned  difficulty  most 
frequently  arises.  Consider  the  type  of  tip  shown  in  figure  4a,  whose  plan- 
form  is  a  semicircle.  The  trailing  edge  is  well-defined  by  a  sharp  corner 
out  to  the  beginning  of  the  tip.  On  the  tip  itself,  the  downstream  side  of 
a  "section"  curve  has  a  finite  radius  of  curvature  which  approaches  zero  at  ( 

the  point  A.  Should  the  trailing  edge  end  at  A  or  should  it  continue  over 
the  tip  region  despite  the  fact  that  there  is  no  sharp  corner?  If  the  i 

"section"  curves  on  the  tip  region  had  sharp  corners,  presumably  the  trailing 
edge  would  continue  into  the  tip  region  all  the  way  to  the  point  B.  For  ( 

highly  yawed  flow,  the  point  B  appears  to  be  part  of  the  leading  edge. 

Where  should  the  trailing  edge  end  in  that  case?  The  tip  in  figure  4b  is  a 
half-body  of  revolution  formed  by  rotating  the  symmetric  section  curve  at 
AA1  about  its  symmetry  line.  In  this  case,  ending  the  trailing  edge  at  the 
point  A  would  probably  be  the  choice  of  most  users.  However,  the  tips  in 
figures  4a  and  4b  differ  mainly  in  their  values  of  the  ratio  of  "spanwise" 
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Figure  4.  Wing  planforms  showing  various  tip  geometries. 


extent  to  "streamwise"  extent.  For  the  "squared-off"  tip  shown  In  figure  4c 
agreement  to  terminate  the  trailing  edge  at  the  point  A  would  be  virtually 
unanimous.  Nevertheless,  the  question  arises  as  to  what  exactly  does  happen 
on  the  tip  Itself.  This  type  of  tip  occurs,  for  example,  at  the  edge  of 
deflected  flaps.  Objections  of  the  sort  mentioned  here  are  basic  to  the 
potential -flow  model  and  do  not  depend  on  the  particular  Implementation  used 
to  produce  an  actual  program. 

One  "answer"  to  the  above  is  that  certain  viscous  effects  are  important 
at  wing  tips,  and  potential  flow  Is  not  expected  to  apply  in  that  region.  The 
"tip  vortex"  leaves  the  wing  well  forward  of  the  trailing  edge  with  a  finite 
diameter  (see  Appendix  B  )  In  contradiction  to  the  potential  flow  model.  Thus, 
the  assumed  potential  flow  model  treats  wing  tips  in  an  approximate  fashion 
and  is  not  applicable  to  very  low  "aspect  ratios". 

A  wing-fuselage  junction  (figure  5a)  Is  another  Important  application 
where  the  trailing  edge  must  end  at  point  A.  It  would  make  little  sense  to 
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Figure  5.  Examples  of  terminating  trailing  edges,  (a)  Wing-fuselage 
intersection,  (b)  A  tip  tank. 


continue  the  trailing  edge  downstream  along,  say,  the  line  AB.  However, 
the  trailing  vortex  wake  Intersects  the  fuselage  along  AB  and  must  do  so 
without  numerical  problems.  The  question  arises  of  what  happens  to  the 
"bound  vortlclty"  at  a  wing-fuselage  junction,  but  that  is  as  much  a  problem 
of  Implementation  as  a  problem  In  the  basic  formulation  (see  section  6.8). 

A  situation  with  elements  of  both  the  above  is  a  wing  with  a  tip  tank 
(figure  5b).  Depending  on  Its  size,  the  tank  may  be  considered  a  small 
fuselage  or  a  big  wing  tip.  Unlike  the  usual  situation  for  a  fuselage,  the 
flow  about  the  tip  tank  has  no  right-and-left  symmetry,  and  there  is  vorticity 
trailing  downstream  from  the  tip  tank,  which  must  be  accounted  for. 

There  are  certainly  other  situations  where  the  details  of  the  potential- 
flow  model  of  three-dimensional  lift  are  unclear.  The  examples  of  this 
section  simply  serve  to  Illustrate  that  such  basic  problems  exist,  regardless 
of  the  particular  Implementation  used  to  reduce  the  model  to  practice.  The 
Implementations  of  course  lead  to  problems  of  their  own. 
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•:  .0  GENERAL  FEATURES  OF  THl  NLTHOD  OF  SOiifTfON 
6.1  The  Method  for  Nonllftlng  Three-Dimensional 

References  1  and  2  review  the  long-tern  effort  of  ».hr  author  and  his 
colleagues  In  the  field  of  potential -flow  circulation.  Ainono  the  methods 
described  are  those  for  lifting  two-dlmens lonal  flows,  and  noolfftlrn  c ?»•.*«"’- 
dimensional  flows.  The  latter  Is  described  In  swrewhat  greater  dataO  «t 
reference  3.  This  nonllftlng  method  forms  tae  Sa:.  :  an  which  Is  boll*:  the 
lifting  method  to  be  described  here.  By  way  of  introduction,  the  nonllf^ng 
program  Is  outlined  briefly  here,  buc  the  references  t*re  relied  or.  t>:.  supply 
all  details. 

All  the  potential  flow  methods  of  references  1  and  Z  are  bused  or.  «t  dis¬ 
tribution  of  source  density  over  the  surface  of  the  body  ab.vjt  v.hith  flew  Is 
to  be  computed.  The  normal  component  of  fluid  velocity  If.  given  on  the 
surface  of  the  body.  Usually  the  normal  velocity  fs  ?ero.  Application  of  t.hP 
norma  I -velocity  boundary  condition  yields  an  Integral  equation  for  the  distri¬ 
bution  function  of  the  source  density,  where  the  domain  of  Integration  Is  the 
body  surface.  Once  this  equation  Is  solved  for  the  source  distribution,  flow 
velocities  both  on  and  off  the  body  surface  can  ne  calculated.  Implement Inn 
this  method  for  the  computer  requires  an  approximate  representation  -;f  tho 
body  surface  and  a  numerical  Integration  routine. 

In  the  nonllftlng  program  of  reference  3,  the  body  Is  specified  to 
computer  by  a  set  of  points,  which  presumably  He  exactly  cr  the  body  surface. 
These  points  are  associated  Into  groups  of  f.nir  "adjacent"  fno4nts  and  a  least- 
squares  plane  passed  through  them.  The  four  pcinti.  are  Jher.  projected  Into 
this  plane  to  form  the  corners  of  a  plane  quadrilateral  surface  element.  Whe* 
this  process  Is  completed  for  all  of  the  points,  the  nody  surface  Is  approx- 
mated  by  a  set  of  plane  quadrilaterals.  A  h/pothetfeol  example  is  snowr*  in 
figure  6.  Because  of  the  process  of  projection,  the  edges  of  adjacent  elements 
may  be  not  quite  coincident,  but  errors  from  this  source  are  small  compared  to 
errors  from  the  other  numerleil  approximations  Inherent  In  the  method. 
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w features  of  the  method  of  approximating  the  body  surface  fr&  of 
importance  to  the  lifting  application.  The  points  defining  the  body  are  Input. 
In  such  an  order  that  they  define  a  family  of  approximately  parallel  nerves 
jyiry  in  the  body  surface.  These  curves,  which  hav*  some  of  the  features  of 
surface  coordinates  have  been  designated  "N-lines,”  as  shown  1r.  figure  6. 

{In  reference  3  the  designation  "colirnm"  Is  used  Instead  of  '’N-llne."  Oath 
have  the  same  meaning.)  First  all  points  along  a  certain  N-linc  are  input  *r. 
order  from  bottom  to  top,  and  then  the  same  is  done  for  the  adjacent  JV-Iine  to 
the  right.  Two  adjacent  N-lincs  bound  a  "strip’  of  elements  of  approximately 
i.r/.istant  width.  The  elements  are  general  quadrilaterals  and  do  not  necessarily 
nave  two  parallel  sides  or  two  r.idas  of  equal  length.  As  6  logical  device  a 
.iiawber  of  N-lines  can  be  associated  Into  a  "section."  Often  a  section  Is 
simply  an  enl<re  body,  but  separate  sections  are  often  used  tc  represent 
ieometrically  different  pe^ts  of  the  same  body;  for  example.  *  wing  and  a 
fuselage.  £lso  sections  are  used  to  concentrate  elements  ir  certain  /eglons 
Of  ...  oo dy.  Logically,  the  concept  of  a  section  means  only  that  the  last  (or 
first)  ‘••line  of  the  section  <s  not  associated  with  the  next  (or  previous) 
H-llw  *.o  form  ?.  strip  of  elements. 


s urf see  e I ewer.vs . 
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On  each  element  one  point  is  selected  where  the  normal  velocity  boundary 
condition  is  to  be  applied  and  where  flow  velocities  are  to  be  computed.  This 
point,  which  Is  designated  the  control  point  of  the  element,  has  been  defined 
various  ways  In  the  past  but  currently  Is  Identified  with  the  centroid  of  the 
element.  Formulas  have  been  derived  that  give  the  components  of  velocity 
Induced  at  a  general  point  In  space  by  a  unit  value  of  source  density  on  a 
general  quadrilateral  element.  These  formulas  allow  the  velocities  induced 
by  the  elements  on  each  other's  control  points  to  be  calculated.  Equating  the 
normal  velocity  Induced  by  all  elements  at  each  control  point  to  the  negative 
of  the  normal  component  of  the  onset  flow  (for  the  case  of  zero  total  normal 
velocity)  yields  a  set  of  linear  algebraic  equations  for  the  values  of  source 
density  on  the  elements.  Once  these  are  solved,  flow  velocities  can  be 
computed  at  the  centroids  and  at  any  selected  point  in  the  flow  field.  For 
the  lifting  application  it  Is  important  to  point  out  that  the  onset  flow  need 
not  be  a  uniform  stream.  Moreover,  solutions  for  several  onsets  may  be 
obtained  simultaneously.  The  onset  flow  affects  only  the  right  side  of  the 
linear  eauations  for  the  source  density  not  the  coefficient  matrix.  Thus,  if 
a  direct  matrix  solution  Is  employed,  several  onset  flows  may  be  treated  in 
nearly  the  same  computing  time  as  a  single  onset  flow. 

6.2  Surface  Elements  for  the  Lifting  Case 

A  lifting  body  and  Its  trailing  vortex  wake  are  approximated  by  quadri¬ 
lateral  surface  elements  in  a  manner  very  similar  to  that  described  in 
reference  3  for  a  nonllftlng  body.  The  approximation  procedure  is  outlined 
here  with  emphasis  on  the  differences  from  the  nonlifting  case. 

As  pointed  out  In  section  5.3,  certain  portions  of  a  general  aerodynamic 
configuration  do  not  have  well-defined  trailing  edges  and  are  not  normally 
thought  of  as  having  their  own  bound  vorticlty;  e.g.,  a  fuselage.  These 
portions  are  denoted  nonlifting  portions  to  signify  that  they  do  not  possess 
Independent  bound  vorticlty  and  that  a  Kutta  condition  is  not  applied  on  them. 
However,  In  general,  the  fluid  exerts  nonzero  pressure  forces  on  nonlifting 
portions  due  to  interference  pressures  from  other  nearby  portions  of  the 
configuration  and  due  to  extentlonc  of  the  bound  vorticlty  from  lifting  portions 
(see  section  6.8).  Nonlifting  portions  are  approximated  by  general  plane 
quadrilateral  elements  In  exactly  the  same  way  as  in  the  nonlifting  method  of 
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reference  3.  In  the  main  calculation  such  elements  have  source  density  but 
not  vortlclty..  The  organization  of  the  Input  data  by  sections  (see  above)  is 
a  natural  way  of  isolating  lifting  and  nonlifting  portions. 

Portions  of  a  general  configuration  that  possess  definite  trailing  edges 
(usually  sharp  corners)  and  contain  bound  vortlclty  are  denoted  lifting  portions. 
The  most  frequently  occurring  application  with  both  lifting  and  nonllfting 
portions  Is  a  wing-fuselage.  Accordingly,  this  configuration  Is  used  as  an 
Illustrative  example  in  figure  7.  On  a  lifting  portion  the  N-llnes  are 
approximately  in  the  freestream  direction.  On  each  N-llne  points  are  Input 
beginning  at  the  trailing  edge,  continuing  around  a  "section  curve"  of  the 
wing,  returning  to  the  trailing  edge,  and  proceeding  downstream  to  define  the 
trailing  vortex  wake.  The  wake  may  be  defined  as  far  downstream  as  desired. 
Provision  has  been  made  to  consider  the  last  element  of  the  wake  semiinfinite 
so  that  wake  definition  may  be  terminated  at  any  point  aft  of  which  the  wake 
curvature  in  the  stream  direction  may  be  neglected.  Usually  a  lifting  portion 
such  as  a  wing  Is  considered  a  single  lifting  section,  but  It  may  be  divided 


Figure  7.  Typical  lifting  configuration. 
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into  several  lifting  sections  if  desired. Within  each  lifting  section  all  N-lines 
must  contain  the  same  number  of  Input  points.  Points  on  adjacent  N-lines  of  a 
lifting  section  are  associated  to  form  surface  elements.  The  set  of  elements 
formed  from  points  on  a  pair  of  adjacent  N-lines  is  denoted  a  "lifting  strip" 
of  elements.  The  strip  contains  elements  both  on  the  body  and  in  the  wake. 
Although  two  adjacent  N-lines  are  not  quite  parallel  in  general,  they  are 
nearly  parallel  in  most  cases. 

Elements  of  lifting  sections  are  taken  as  plane  trapezoids.  Each  of  the 
two  parallel  sides  is  formed  from  two  input  points  on  the  same  N-line.  Thus 
the  parallel  sides  are  approximately  along  the  N-lines.  Of  course,  in  the 
general  case  the  four  input  points  that  are  associated  to  form  an  element  do 
not  even  lie  in  the  same  plane,  much  less  form  a  trapezoid.  They  must  be 
"adjusted"  to  do  this.  In  the  nonlifting  program  of  reference  3  the  input 
points  are  adjusted  to  lie  in  the  same  plane  but  not  to  be  trapezoidal.  Thus, 
the  "adjustment"  required  is  somewhat  more  for  lifting  elements  than  for  non- 
lifting.  Adjacent  elements  have  two  Input  points  in  common,  but  the  adjustment 
that  these  points  are  subject  to  is  usually  different  for  the  two  elements. 

Thus,  In  general,  after  adjustment  the  sides  of  adjacent  elements  are  not 
coincident,  and  there  are  gaps  between  the  elements.  Such  gaps  exist  for  both 
lifting  and  nonlifting  elements.  For  the  nonlifting  case  the  unimportance  of 
the  gaps  is  discussed  in  references  1  and  3.  For  lifting  elements  the  gaps  are 
presumably  greater  than  for  nonlifting  elements,  but  it  seems  that  in  both  cases 
the  gaps  should  have  the  same  order  of  magnitude.  Thus,  errors  from  this  source 
should  be  unimportant.  It  is  pointed  out  in  references  1  and  3  that  for  some 
bodies  the  gaps  between  elements  vanish.  For  lifting  bodies  the  1.v>or co : ■  t  case 
for  which  this  occurs  is  an  untwisted  wing,  possibly  swept  and  tapored,  I  /ing 
the  same  airfoil  section  at  all  spanwise  locations. 

The  centroids  of  the  elements  are  used  as  control  points.  Thus,  for  each 
lifting  strip  the  locus  of  control  points  is  approximately  midway  between  the 
two  N-lines  used  to  generate  the  strip.  Elements  of  lifting  strips  have 
source  densities  whose  strengths  are  determined  to  give  zero  (or  prescribed) 
normal  velocity  at  the  control  points. 
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6.3  Bound  and  Trailing  Vorticity 


In  addition  to  the  source  densities  on  the  elements,  lifting  portions 
also  possess  a  distribution  of  bound  vorticity.  As  pointed  out  In  section 
5.2,  the  form  of  the  bound  vorticity  uniquely  determines  the  strength  distri¬ 
bution  of  the  trailing  vorticity,  which  lies  along  the  input,  wake.  The  form 
assumed  for  the  bound  vorticity  contains  a  number  of  adjustable  parameters 
equal  to  the  number  of  lifting  strips  on  that  lifting  portion.  The  values 
of  these  parameters  are  determined  by  applying  a  Kutta  condition  at  the 
trailing  edge  segment  (figure  7)  of  each  lifting  strip.  The  simplest  form 
of  the  bound  vorticity  distribution  utilizes  a  set  of  individual  distribu¬ 
tions,  each  of  which  is  nonzero  only  on  one  lifting  strip.  The  complete 
distribution  consists  of  a  linear  combination  of  these  Individual  distribu¬ 
tions,  each  of  which  is  nonzero  on  a  different  lifting  strip.  The  combination 
constants  of  the  linear  combination  are  the  required  adjustable  parameters. 
This  is  the  type  of  distribution  used  In  the  present  method.  Other  existing 
methods  (references  4,  5,  6,  and  7)  also  use  this  type  of  distribution. 

The  value  of  the  parameter  multiplying  the  distribution  associated  with  a 
particular  lifting  strip  represents  the  strength  of  the  bound  vorticity  at 
the  "spanwise"  location  of  that  strip.  Thus,  as  expected,  the  "spanwise" 
variation  of  bound  vorticity  is  determined  by  the  Kutta  condition.  More 
precisely  the  "spanwise"  variation  cf  vorticity  from  one  lifting  strip  to 
another  is  determined  by  the  Kutta  condition.  The  "spanwise"  variation  of 
vorticity  within  the  small  but  finite  span  of  each  individual  lifting  strip 
Is  basically  a  question  of  the  order  of  accuracy  of  a  numerical  integration 
(see  below  for  the  options  of  the  present  method). 

Even  If  the  bound  vorticity  is  of  the  type  mentioned  above,  various 
forms  of  this  vorticity  are  possible.  In  addition,  the  "chordwise"  or 
"streamwise"  variation  of  vorticity  on  a  "section  curve"  at  a  particular 
"spanwise"  location  may  be  chosen  at  will.  In  the  limit  where  an  infinite 
number  of  surface  elements  are  used  to  approximate  the  body.  It  appears  that 
the  calculated  flow  velocities  are  independent  of  the  assumptions  made  con¬ 
cerning  bound  vorticity.  However,  for  practical  element  numbers,  the  form 
assumed  for  the  bound  vorticity  and  its  "chot ,  'ise"  variation  have  an 
appreciable  effect  on  the  accuracy  of  the  sol;  -■  i.  The  methods  of 
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references  4,  5,  6,  and  7  all  use  the  same  form  for  the  bound  vorticlty, 
which  consists  of  concentrated  vortex  filaments  lying  In  the  camber  surface 
of  the  wing.  Some  details  are  illustrated  in  figure  8a,  which  shows  a 
single  N-llne  representing  a  section  curve  of  the  wing.  An  equal  number  of 
elements  is  placed  on  the  upper  and  lower  surfaces.  The  input  points  defining 
the  elements  are  arranged  so  that  a  pair  of  points,  one  on  the  upper  surface 
and  one  on  the  lower,  lie  nearly  on  the  same  perpendicular  to  the  mean 
camber  surface.  The  bound  vorticlty  filaments,  which  appear  as  points  in 
figure  8a,  lie  midway  between  corresponding  points  on  the  upper  and  lower 
surface.  This  arrangement  maximizes  the  distance  of  the  vortex  filaments 
from  the  wing  surface  and  presumably  reduces  numerical  problems  associated 
with  the  flow  singularities  at  the  filaments.  Thus,  in  general  the  number 
of  vortex  filaments  is  one  less  than  half  the  number  of  surface  elements  in 
the  lifting  strip,  although  in  certain  formulations  some  vortices  may  be 
given  zero  strength.  The  strengths  of  the  bound  vortex  filaments  are  main¬ 
tained  constant  over  the  "span"  of  each  individual  lifting  strip.  Thus, 


Figure  8.  Representation  of  the  bound  vorticity  by  concentrated  vortex 
filaments  lying  in  the  mean  camber  surface,  (a)  A  section 
curve  of  the  wing,  (b)  The  complete  three-dimens  tonal 
vortex  pattern. 
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the  trailing  vorticity  is  also  in  concentrated  filaments.  Forward  of  the 
trailing  edge  these  lie  In  the  mean  camber  surface  beneath  the  edges  of  the 
strip,  i.e.,  midway  between  the  portions  of  the  N-lines  on  the  upper  and 
lower  surfaces  of  the  wing.  Downstream  of  the  trailing  edge  the  trailing 
vortex  filaments  He  along  the  N-llnes  defining  the  assumed  wake.  A  view  of 
the  entire  three-dimensional  arrangement  is  shown  in  figure  8b.  The  formula¬ 
tions  of  the  references  use  different  "chordwise"  variations  of  the  vortex 
strengths.  Reference  4  presents  results  for  a  distribution  of  zero  strength 
from  0%  to  20%  chord  and  from  80%  to  100%  chord.  From  20%  to  80%  chord  the 
distribution  is  constant.  However,  both  reference  4  and  the  subsequent 
development  of  the  method  presented  in  reference  5  recommend  use  of  a 
"chordwise"  vorticity  variation  approximately  the  same  as  the  "chordwise" 
lift  distribution.  In  a  practical  case  this  last  might  be  determined  from 
linear  theory  or  might  be  estimated  from  results  for  similar  wings.  Quite 
different  are  the  distributions  used  in  references  6  and  7.  Apparently, 
reference  6  uses  a  vortex  strength  proportional  to  the  local  thickness  of 
the  airfoil  section,  while  reference  7  uses  a  strength  proportional  to  the 
square  root  of  the  local  thickness.  Since  exact  solutions  are  not  available 
and  experimental  results  are  affected  by  viscosity,  compressibility,  and 
testing  error,  the  results  of  these  calculations  must  be  judged  largely  on 
their  "reasonableness,"  e.g.,  lack  of  extraneous  wiggles,  etc. 

The  present  method  uses  a  completely  different  form  for  the  bound 
vorticity.  Instead  of  concentrated  vortex  filaments  interior  to  the  wing, 
there  is  a  finite-strength  sheet  of  vorticity  on  the  surface  of  the  wing, 
i.e.,  the  vorticity  lies  on  the  quadrilateral  surface  elements.  The  nature 
of  the  singularity  is  thus  reduced  from  .  ^ne  singularity  to  a  surface 
singularity.  Some  features  of  this  formul,  :  -p  are  Illustrated  in  figure  9 
which  may  be  compared  with  figure  8.  The  "chordwise"  variation  of  the  surface 
vorticity  strength  may  be  chosen  at  will.  In  the  present  method  the  strength 
is  taken  as  constant  all  around  the  airfoil  section.  This  choice  was  influ¬ 
enced  by  requirements  of  simplicity  and  by  the  fact  that  constant-strength 
surface  vorticity  gives  good  results  In  two-dimensional  cases  (see  below). 

The  variation  of  vorticity  over  the  "span"  of  a  lifting  strip  of  elements  has 
two  options:  constant  and  linear.  In  the  former  option  the  "spanwise"  vari¬ 
ation  of  vorticity  over  the  wing  is  a  step  function  (figure  10a)  whose  values 
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SURFACE 


Figure  9.  Representation  of  the  bound  vorticity  by  a  finite-strength 

vortlclty  distribution  lying  on  the  wing  surface,  (a)  A  section 
curve  of  the  wing,  (b)  The  complete  three-dimensional  vorticity 
pattern  using  a  step  function  spanwise  variation,  (c)  The 
complete  three-dimensional  vorticity  pattern  using  a  piecewise 
linear  spanwise  variation. 
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Figure  10.  Two  forms  of  the  spanwise  variation  of  bound  surface 
vorticity.  (a)  Step  function,  (b)  Piecewise  linear. 
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are  determined  by  the  Kutta  condition.  This  form  of  the  bound  vorticity 
has  the  advantage  of  simplicity  and  does  not  require  special  handling  at 
the  end  of  a  lifting  section,  e.g.,  a  wing  tip.  However,  the  trailing 
vorticity  takes  the  form  of  concentrated  vortex  filaments  along  the  N-lines 
(figure  9b).  This  situation  can  be  avoided  by  using  a  linear  vorticity 
variation  over  the  span  of  the  lifting  strip.  In  this  case  the  trailing 
vorticity  takes  the  form  of  a  vortex  sheet  over  the  surface  of  the  strip, 
i.e.,  over  the  surface  elements  (figure  9c).  If  the  vorticity  distribution 
were  exactly  continuous  at  the  edges  of  the  strips,  i.e.,  at  the  N-lines, 
there  would  be  no  vortex  filaments  on  the  N-lines.  This  is  not  possible  in 
general  because,  as  mentioned  in  section  6.2,  the  edges  of  adjacent  elements 
are  not  quite  coincident.  Thus,  there  are  small  geometrical  discontinuities 
in  the  vortex  sheet  along  the  N-lines.  It  is  thus  not  worthwhile  to  attempt 
to  determine  the  "spanwise"  rate  of  change  of  vorticity  over  a  strip  from  a 
condition  of  continuity  of  strength  along  the  N-lines.  Moreover,  this  type 
of  variation  leads  to  serious  numerical  difficulties  (reference  8).  Instead 
the  spanwise  rate  of  change  on  a  strip  is  determined  from  a  centered  parabolic 
fit  over  values  of  bound  vorticity  at  the  midspan  of  three  consecutive  strips 
and  strict  continuity  of  strength  at  the  N-lines  is  obtained  only  if  the 
"spanwise"  variation  is  truly  parabolic.  However,  the  discontinuity  is  of 
high  order,  and  the  vortex  sheet  may  be  considered  continuous  to  within  the 
order  of  the  overall  approximation.  In  this  option  the  "spanwise"  variation 
of  vorticity  Is  a  piecewise  linear  function  as  shown  in  figure  10b.  The 
trailing  vorticity  continues  as  a  sheet  into  the  wake,  so  that  the  velocity 
has  the  desired  behavior  of  discontinuity  across  the  wake.  The  behavior  does 
not  occur  If  the  wake  is  composed  of  concentrated  filaments  as  it  is  in  the 
methods  of  the  references  and  In  the  above  "step  function"  option  of  the 
present  method.  The  chief  disadvantage  of  the  "piecewise  linear"  option  is 
that  special  handling  is  required  at  the  first  and  last  lifting  strips  of  a 
section  to  determine  the  "spanwise"  rate  of  change  of  vorticity  (section  7.11). 
Moreover,  in  most  cases  that  have  been  run  with  the  present  method  using 
both  options  for  the  bound  vorticity,  the  calculated  results  are  not  very 
different. 

The  accuracy  to  be  obtained  using  various  forms  for  the  bound  vorticity 
may  be  investigated  by  considering  the  two-dimensional  case  for  which  exact 
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analytic  solutions  are  available.  Indeed  this  is  a  very  natural  procedure 
because  the  essential  three-dimensional  feature  is  the  "spanwise"  variation 
of  vortlcity  which  is  determined  by  the  Kutta  condition.  The  form  of  the 
bound  vortlcity  and  its  assumed  "chordwlse"  variation  have  direct  two- 
dimensional  analogies,  which  are  very  similar  numerically  to  what  is  being 
calculated  in  three  dimensions.  The  two-dimensional  cases  are  obtained  by 
simply  considering  the  "section  curves"  of  figures  8a  and  9a  as  two- 
dimensional  airfoils.  The  cases  were  run  with  the  rather  small  element 
numbers  that  are  characteristic  of  the  three-dimensional  case  rather  than 
the  much  larger  element  numbers  that  are  available  in  two  dimensions  to 
obtain  very  high  accuracy.  Two  cases  are  presented  here  that  illustrate 
different  aspects  of  the  situation. 

The  first  case  is  a  Karman-Trefftz  airfoil,  for  which  coordinates  of 
points  on  the  body  may  be  obtained  very  accurately  using  analytic  expressions. 
A  rather  extreme  geometry  was  chosen  so  that  differences  in  the  solutions 
could  be  seen  more  easily.  The  airfoil  is  8.2  percent  thick,  has  a  9° 
trail inq-edge  angle  and  the  rather  large  camber  value  of  24  percent.  A  sketch 
of  the  shape  is  given  In  figure  11.  Calculations  were  performed  for  an  angle 
of  attack  of  1.205°,  The  exact  solution  from  the  well-known  formulas  gives 
a  lift  coefficient  of  3.37.  Using  50  surface  elements,  calculations  were 
performed  with  a  constant-strength  surface  vortlcity,  as  is  done  in  the 
present  method,  and  also  with  interior  vortex  filaments  whose  strength  is 
proportional  to  the  local  airfoil  thickness,  as  Is  done  in  the  method  of 
reference  6.  The  calculated  surface  pressure  distributions  are  compared  with 
the  exact  solution  In  figure  11.  Neither  calculated  result  is  very  good 
because  of  the  extreme  geometry  and  the  limited  element  number.  However,  the 
error  for  the  surface  vorticity  approach  is  about  half  the  error  for  the 
Interior  vortex  filament  approach.  The  "wiggles"  in  the  solution  generated 
from  the  interior  vortex  filaments  are  not  due  to  inaccuracies  in  the  points 
defining  the  airfoil.  These  points  are  exact.  The  "wiggles"  are  apparently 
due  to  changes  In  element  lengths  along  the  surface.  Adjacent  elements  differ 
in  length  by  no  more  than  25  percent,  which  appears  quite  reasonable.  The 
solution  obtained  from  the  surface  vortlcity  does  not  respond  to  this  situa¬ 
tion  and  is  perfectly  smooth. 
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Figure  11.  Surface  pressure  distributions  on  a  Karman-Trefftz  airfoil 
of  large  canter  at  1.205  degrees  angle  of  attack. 


calculated  solutions 

-  EXACT 
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Figure  12.  Surface  pressure  distributions  on  a  conventional  airfoil  section  at  6.9  degrees  angle  of  attack. 


The  second  case  Is  the  conventional  airfoil  section  shown  in  figure  12. 

The  coordinates  of  the  points  defining  this  airfoil  were  obtained  by  procedures 
usual  in  design  applications,  and  the  result  is  that  the  point  distribution  Is 
not  absolutely  smooth  but  contains  small  Irregularities.  Calculations  were 
performed  with  32  surface  elements,  Figure  12  shows  the  points  defining  the 
airfoil  and  the  locations  of  the  15  Interior  vortex  filaments  that  were  used 
in  the  calculations  with  strengths  proportional  to  local  thickness.  Calcula¬ 
tions  were  also  performed  using  the  constant-strength  surface  vorticlty  of 
the  present  method.  Surface  pressure  distributions  calculated  by  the  two 
methods  are  compared  with  a  very  accurate  conformal -mapping  solution  in 
figure  12.  The  surface  vorticlty  approach  Is  unaffected  by  any  Irregularities 
of  the  points  and  its  results  agree  very  well  with  the  accurate  solution. 

In  fact  the  point  distribution  of  figure  12  Is  the  one  used  with  the  present 
method  to  produce  the  three-dimensional  results  of  figure  42.  The  pressure 
distribution  calculated  by  the  approach  based  on  Interior  vortex  filaments 
has  rather  severe  "wiggles"  and  also  has  a  systematic  error  in  pressure 
level  so  that  the  lift  coefficient  obtained  by  integrating  the  pressures 
differs  from  the  exact  value  by  20  percent. 

From  these  two  examples  and  others  that  have  been  run,  it  is  concluded 
that  the  representation  of  the  bound  vorticlty  by  finite-strength  surface 
vorticlty  Is  superior  to  the  representation  by  interior  vortex  filaments. 

The  former  is  far  less  sensitive  to  inaccuracy  of  the  input  data  and  tends 
to  give  a  more  accurate  solution  even  when  the  data  is  smooth. 

6.4  Use  of  a  Dipole  Distribution  to  Represent  Vorticlty 

From  the  previous  section  It  can  be  seen  that  in  the  present  method  the 
bound  and  trailing  vorticity  are  represented  by  a  general  surface  distribution 
of  vorticity,  possibly  with  concentrated  vortex  filaments  at  the  edges. 

Formulas  that  express  the  velocity  Induced  by  such  a  vorticity  distribution 
are  required.  Derivation  of  such  expressions  is  complicated  by  the  fact  that 
the  surface  vorticity  strength  is  a  vector  that  varies  in  both  magnitude  and 
direction.  Furthermore,  care  must  be  taken  to  insure  that  the  vorticity 
distribution  gives  rise  to  a  potential  flow,  i.e.,  that  the  individual 
infinitesimal  vortex  lines  either  form  closed  curves  or  go  to  infinity.  Use 
of  a  surface  dipole  distribution  circumvents  these  complications,  because 
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the  dipole  strength  is  a  scalar  and  any  arbitrary  diw>1e  distributee?*  gives 
rise  to  a  potential  flow.  A  general  result  *; f v i *«g  the  rclatfonsMo  between 
dipole  sheet  and  a  vortex  sheet  is  given  in  Appendix  ft.  It  my  be  sunriA**1;r*»ii 
as  follows:  A  variable-strength  dipole  sheet  is  equivalent  to  the  swc  tJ : 

(1)  a  variable-strength  vortex  sheet  on  the  sime  svrlace  as  the  dioole  sheet 
whose  vorticity  has  a  direction  at  right  angles  to  th.e  gradient  of  the  Jipoie 
strength  and  a  magnitude  equal  to  the  magnitude  of  this  g**adi**nt,  and  (2)  a 
concentrated  vortex  filament  around  the  edge  of  the  shpet  whose  strength  U 
everywhere  equal  to  the  local  edge  value  of  ‘pole  strength,  fhis  relation., 
which  is  a  straightforward  generalization  of  tie  well-know*  two-dlmeoslr,*  a', 
result,  does  not  appear  explicitly  in  the  literatus*.  Its  plausibility  was 
discussed  early  in  the  present  work  by  the  a  ithor,  A.M.O.  smith,  and  P.8.5. 
Lissaman.  The  proof  of  this  relation  in  Appendix  A,  which  wa-  originally 
outlined  by  the  author  in  reference  9,  is  apparently  the  first.  A  late* 
derivation  is  contained  in  reference  10.  In  the  preser.x  «?iethou  a^l  formulas 
are  derived  in  terms  of  dipole  distributions  and  thfr  above  relationship  Is  used 
to  interpret  this  situation  In  terms  of  the  more  physically  significant  vorticity. 
In  particular  "chordwlse"  dipole  variation  is  equivalent  to  "snamrise"  vorticity 
and  "spanwise"  dipole  variation  to  "chordwlse’’  vorticity.  Also,  if  a  dipole 
sheet  terminates  with  a  nonzero  strength,  it  results  in  a  concentrated  vortex 
filament. 


6.5  The  Kutta  Condition 

It  is  an  Interesting  and  Important  fact  that  the  "physical*  Kutti  condition 
of  finite  velocity  at  the  trailing  edge  cannot  be  applied  In  a  general  numerical 
procedure  for  calculating  flow.  This  Is  true  In  both  two  dimensions  and  three 
dimensions.  If  the  general  solution  could  be  written  down  in  explicit  analytic 
form,  as  is  possible  in  a  few  simple  two-dimensional  cases,  then  the  appropriate 
parameters  could  be  adjusted  to  eliminate  the  singular  terms  in  the  expression 
for  surface  velocity.  However,  in  a  numerical  solution  there  is  no  true  singu¬ 
larity,  and  a  condition  of  finiteness  without  specifying  a  definite  value  cannot 
determine  specific  values  of  a  parameter.  Accordingly,  the  kutta  condition  is 
applied  by  Indirect  means.  What  Is  done  Is  to  deduce  another  property  of  the 
flow  at  the  trailing  edge  that  Is  a  direct  consequer.ee  of  the  finiteness  of  vel¬ 
ocity  and  to  use  this  related  property  as  "the  Kutta  condition."  Various  properties 
may  be  derived.  Some  are  strictly  valid  onl>  fer  the  true  flow  (limit  of  infinite 
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wbcr  o»  elements)  *nd  are  applies  u>  a  case  of  finite*  element  number  as  an 
approxi  ..cfoo.  C-th nrz  napper.  to  bo  t'-ue  for  finite  element  number,  and  still 
others  hav°  dl»ferer*t  forsf.i  ;r*  cas^s  of  infinite  and  of  finite  element  numbers. 
In  general,  conditions  r.jnmot  be  exactly  at  the  trailing  edge  If  a 

finite  number  of  elements  1*-  used  (except  In  the  sense  that  quantities  can 
be  extrapolated  to  the  trailing  *doe).  Thus,  "the  Kutta  condition"  is  applied 
a  small  distance  away  fran  the  trailing  edge,  and  determining  an  appropriate 
value  for  this  distance  arid  its  effect  or  the  solution  is  part  of  the  problem. 
The  situation  can  Le  effected  by  t fie  fact  that  some  flow  conditions  at  the 
trading  edge  are  extrersely  local,  and  their  values  are  quite  different  even 
a  small  distance  away.  Such  very  local  conditions  cannot  be  applied  to 
case:*  of  reasonable  element  numbers. 


Some  related  properties  that  nay  be  deduced  from  tha  Kutta  condition 

arp  as  follows: 

Two-Pimens ions: 

(a)  A  streamline  of  the  flow  leaves  the  trailing  edge  with  a  direction 
along  the  hi sector  uf  the  tralling-edge  angle. 

•b)  As  toe  trailing  edge  is  approached  the  surface  pressures  (velocity 
magnitudes)  on  the  upper  and  lower  surfaces  have  a  common  limit, 
which  equals  stagnation  pressure  (zero  velocity)  if  the  trailing- 
edqe  angle  is  oonzero. 

(c)  Tha  source  density  at  the  trailing  edge  is  zero. 


Three-Dimeosions: 

(a)  A  stream  surface  of  the  flow  leaves  the  trailing  edge  with  a 
direction  that  is  known,  or  at  least  can  be  approximated  (see 

belcw) . 

(b)  As  the  trailing  edge  is  approached,  the  surface  pressures  (velocity 
magnitudes)  on  the  upper  and  lower  surfaces  have  a  common  limit. 

(c)  The  source  density  at  the  trailing  edge  is  zero. 

The  example  properties  above  can  be  used  to  apply  the  Kutta  condition  in 
cases  of  finite  element  number.  Property  (a)  In  either  dimensionality 
differs  f»om  the  others  in  that  it  must  be  applied  off  the  body  surface. 
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Points  downstream  of  the  trailing  edge  are  selected  to  be  on  the  stream 
surface  or  streamline  and  directions  normal  to  the  stream  surface  or  stream¬ 
line  are  prescribed.  Then  a  flow  tangency  condition  of  zero  normal  velocity 
is  applied  at  these  points  just  as  If  they  were  control  points  of  rface 
elements.  Selection  of  distances  from  the  trailing  edge  at  which  to  apply 
the  flow  tangency  condition  Is  part  of  the  problem.  Properties  (b)  and  (c) 
are  applied  on  the  body  surface.  Since  the  flow  on  the  body  has  meaning 
only  at  the  control  points,  these  conditions  are  applied  to  flow  quantities 
at  the  control  points  of  the  elements  adjacent  to  the  trailing  edge  on  the 
upper  and  lower  surfaces.  In  two  dimensions  there  are  just  two  such  elements, 
while  in  three  dimensions  there  are  two  elements  on  each  lifting  strip.  It 
might  be  supposed  that  property  (c)  is  applied  by  requiring  source  densities 
on  elements  adjacent  to  the  trailing  edge  to  be  zero.  This  amounts  to  two 
conditions  per  lifting  strip  and  thus  overdetermines  the  problem.  The  best 
that  can  be  done  Is  to  require  that  for  each  lifting  strip  the  values  of 
source  density  on  the  two  elements  adjacent  to  the  trailing  edge  be  equal  in 
magnitude  and  of  opposite  sign.  Similarly,  condition  (b)  is  applied  by 
requiring  that  for  each  lifting  strip  the  magnitudes  of  the  velocity  at 
the  control  points  of  the  two  elements  adjacent  to  the  trailing  edge  be 
equal.  This  Is  done  even  In  two  dimensions  where  the  theoretical  velocity 
of  zero  Is  so  local  that  the  velocity  is  an  appreciable  fraction  of  freestream 
velocity  at  the  control  point  adjacent  to  the  trailing  edge. 

In  applications,  property  (c)  has  not  been  used.  The  methods  of 
references  4,  5,  6  and  7  use  property  (a).  The  present  method  has  the  option 
of  using  either  property  (a)  or  property  (b)  as  "the  Kutta  condition."  If 
property  (a)  Is  used  the  points  where  it  is  to  be  applied  and  the  normal 
vectors  at.  these  points  must  be  furnished  to  the  program  as  input.  Flow 
velocities  are  computed  at  all  control  points  due  to  the  bound  vorticlty 
distribution  associated  with  each  lifting  strip.  Each  of  these  flows  is 
considered  as  an  onset  flow  to  the  body.  Let  the  total  number  of  quadri¬ 
lateral  source  elements  be  N  and  the  number  of  lifting  strips  be  L.  Then 
there  are  L  vorticlty  onset  flows,  each  of  which  consists  of  velocity  com¬ 
ponents  at:  the  N  control  points,  the  L  points  where  property  (a)  is  to 
be  applied  (If  that  option  is  used),  and  any  other  off-body  point  where  flow 
is  to  be  computed.  For  each  onset  flow  a  set  of  N  values  of  source  density 
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on  the  elements  is  obtained  that  gives  zero  normal  velocities  at  the  N 
control  points.  The  same  is  done  for  the  uniform  onset  flow  that  represents 
the  freestream.  As  described  In  section  6.1,  the  values  of  source  density 
are  obtained  as  solutions  of  a  set  of  linear  algebraic  equations  whose 
N  x  N  coefficient  matrix  is  the  same  for  all  L  +  1  onset  flows.  The  onset 
flows  simply  yield  L  +  1  right  sides  for  the  equations.  Using  a  direct 
matrix  solution  all  L  +  1  sets  of  source  density  are  obtained  simultaneously. 
The  desired  source  density  distribution  is  a  linear  combination  of  these 
individual  distributions.  The  constants  in  this  linear  combination  are  the 
L  values  of  bound  vorticlty  associated  with  the  various  lifting  strips,  and 
these  are  determined  from  the  Kutta  condition.  (The  solution  corresponding 
to  the  uniform  stream  enters  with  unit  coefficient.)  Flow  velocities  for 
the  individual  solutions  are  computed  only  for  the  points  used  to  apply  the 
Kutta  condition  —  either  the  control  points  of  the  elements  adjacent  to  the 
trailing  edge  If  property  (b)  Is  used,  or  the  additional  input  points  down¬ 
stream  of  the  trailing  edge  if  property  (a)  is  used.  The  Kutta  condition 
results  in  L  simultaneous  equations  whose  solution  yields  the  desired  L 
values  of  bound  vorticlty.  In  typical  cases  the  number  of  lifting  ?frips  L 
is  10  to  30,  as  contrasted  with  the  number  of  surface  elements  N,  which  is 
300  to  1000.  Thus,  solution  of  the  equations  expressing  the  Kutta  condition 
Is  a  negligible  computation  compared  to  solution  of  the  equations  for  the 
values  of  source  density.  The  values  of  bound  vorticity  are  used  to  compute 
a  single  set  of  N  values  of  source  density  -  the  "combined"  values  -  that 
are  used  to  compute  velocities  at  the  control  points  of  the  elements. 

An  alternative  numerical  procedure  for  implementing  the  application  of 
the  Kutta  condition  is  employed  in  references  6  and  7.  As  mentioned  above, 
the  bound  vorticlty  associated  with  each  lifting  strip  induces  a  velocity  at 
each  control  point.  These  may  be  treated  exactly  the  same  as  the  velocities 
induced  by  the  individual  source  quadrilaterals  (section  6.1),  i.e.,  the  I. 
values  of  bound  vorticlty  may  be  treated  as  additional  unknowns  in  the  equa¬ 
tions  expressing  the  normal  velocity  boundary  condition.  This  yields  N 
linear  equations  In  N  +  L  unknowns.  The  Kutta  condition  supplies  the  addi¬ 
tional  L  equations.  If  the  Kutta  condition  is  expressed  as  property  (a), 
as  is  done  in  references  6  and  7,  the  additional  L  equations  are  linear. 
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As  discussed  in  references  1,  2,  and  3,  the  N  x  N  coefficient  matrix 
due  to  the  source  quadrilaterals  has  a  dominant  main  diagonal  and  is  well 
suited  to  numerical  solution  either  by  direct  solution  or  by  iterative  solu¬ 
tion.  The  additional  L  equations  from  the  Kutta  condition  do  not  have 
dominant  diagonal  terms  so  that  the  (N  +  L)  x  (H  +  L)  matrix  used  in  refer¬ 
ences  6  and  7  is  not  well-conditioned.  However,  suitable  partitioning  of 
this  matrix  (the  partitioning  is  different  in  reference  6  from  that  of  ref¬ 
erence  7)  yields  rapidly  convergent  iterative  solutions.  If  the  property  (b) 
is  used  to  express  the  Kutta  condition,  the  additional  L  equations  are 
quadratic  because  they  are  applied  to  a  vector  magnitude.  (In  two  dimensions 
the  surface  velocity  has  only  one  component,  and  the  equations  derived  from 
property  (b)  are  linear.)  This  might  not  be  a  serious  handicap  in  an  iter¬ 
ative  procedure,  but  it  has  never  been  tried. 

The  relative  advantage  or  disadvantage  of  an  iterative  solution,  like 
that  of  references  6  and  7,  compared  to  a  direct  solution,  like  that  of  the 
present  method,  is  primarily  a  matter  of  computing  time.  The  situation  is 
affected  by  the  particular  computer  being  used  and  by  the  accounting 
algorithm  for  multiple-user  machines.  However,  by  fcr  the  most  important 
considerations  are  the  element  number  N  and  the  type  of  body  about  which 
flow  is  to  be  computed.  A  direct  solution  for  a  set  of  linear  equations 
requires  a  computing  time  proportional  to  N3,  and  this  time  is  independent 
of  the  body.  An  iterative  solution  requires  a  computing  time  proportional 
to  the  product  IN2,  where  I  is  the  number  of  Iterations  needed  for 
convergence.  It  Is  clear  then  that  for  sufficiently  large  N,  the  iterative 
solution  Is  quicker  (assuming  that  I  is  Independent  of  N,  which  appears 
to  be  the  case  in  the  present  application).  Similarly,  for  sufficiently  small 
N  the  direct  solution  Is  quicker.  The  "crossover"  value  of  N,  where  the 
two  methods  are  equal  is  directly  proportional  to  I.  For  simple 
bodies,  such  as  wing-fuselages,  I  Is  approximately  15  and  the  crossover 
value  for  N  Is  perhaps  800  for  an  IBM  370-165.  In  any  event,  the  Iterative 
solution  is  clearly  superior  for  N  =  1000,  and  the  direct  solution  Is 
clearly  superior  for  N  =  500.  For  more  complicated  bodies,  and  particularly 
for  situations  Involving  Interior  flows,  I  Is  considerably  larger,  and  thus 
so  Is  the  crossover  value  of  N.  Such  situations  arise,  for  example,  with 
nacelles  (reference  1)  and  with  bodies  in  a  wind  tunnel  (section  9.4).  If 
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the  estimated  computing  times  are  not  too  different,  the  direct  solution  is 
to  be  preferred,  because  the  time  required  is  predictable.  It  appears  that 
the  most  efficient  procedure  is  one  containing  both  direct  and  iterative 
solut  of  the  linear  equations  as  options.  Inclusion  of  an  iterative 
solution  in  the  present  method  is  a  desirable  future  extension. 

In  the  present  method,  application  of  property  (b)  is  straightforward 
and  requires  no  additional  input.  Its  effectiveness  can  be  judged  simply 
by  the  accuracy  of  the  resulting  calculation,  as  discussed  below.  Applica¬ 
tion  of  property  (a)  (either  in  the  present  method  or  in  the  methods  of 
references  6  and  7)  requires  the  answer  to  two  questions:  How  far  from  the 
trailing  edge  should  it  be  applied?  In  what  direction  with  respect  to  the 
trailing  edge  should  the  point  of  application  be  situated?  The  answer  to 
the  second  question  which  will  be  considered  first,  appears  to  be  related 
to  the  direction  by  which  the  stream  surface  leaves  the  trailing  edge  of 
the  wing.  However,  this  last  turns  out  to  be  false  In  many  applications. 

The  behavior  of  the  vortex  wake  in  the  neighborhood  of  the  trailing  edge 
of  a  three-dimensional  lifting  body  has  been  worked  out  from  basic  principles 
in  reference  11  under  the  assumption  of  Inviscid  potential  flow.  The  results 
are  easy  to  state.  The  only  two  quantities  that  affect  the  situation  are 
the  local  section  lift  coefficient  and  the  local  value  of  the  average  component 
of  velocity  along  the  trailing  edge  (averaged  between  upper  and  lower  surfaces). 
Theoretically,  the  magnitudes  of  these  two  quantities  are  not  important  - 
only  their  signs.  Consider  the  usual  case  when  the  local  section  lift  coef¬ 
ficient  is  positive.  Then  reference  11  states  that  if  the  component  of 
velocity  along  the  trailing  edge  is  outboard,  the  vortex  wake  leaves  the 
trailing  edge  tangent  to  the  upper  surface.  If  this  velocity  component  is 
inboard,  the  sheet  leaves  tangent  to  the  lower  surface.  The  situation  is 
Illustrated  in  figure  13.  If  the  local  section  lift  coefficient  is  negative, 
the  situation  is  reversed. 

The  above  results  mean  that  the  way  in  which  the  vortex  wake  leaves  the 
trailing  edge  depends  on  the  final  flow  solution  and  is  thus  not  known  ahead 
of  time.  On  the  face  of  it  this  is  a  problem.  However,  in  many  practical 
cases  it  is  obvious  which  way  the  flow  at  the  trailing  edge  goes.  In  regions 
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Figure  13.  Theoretical  behavior  of  the  vortex  wake  at  the  trailing 
edqe  of  a  wing,  (a)  Outboard  trailing  edge  velocity, 

(b;  Inboard  trailing  edge  velocity. 

where  there  is  some  doubt  the  flow  component  along  the  trailing  edge  is 
probably  small  compared  to  freestream.  This  situation,  which  occurs  rather 
often  In  applications,  has  an  important  effect  on  the  application  of  the  Kutta 
condition. 


Reference  11  proves  that  for  any  outboard  tral ling-edge  velocity,  no 
matter  how  small,  the  wake  is  tangent  to  the  upper  surface  as  shown  in 
figure  13a.  Similarly,  for  any  inboard  velocity,  no  matter  how  small,  the 
wake  Is  tangent  to  the  lower  surface,  as  shown  In  figure  13b.  On  the  other 
hand,  it  is  physically  evident  that  a  small  change  in  conditions  at  the 
trailing  edge  gives  a  small  change  in  the  wake  position  a  finite  distance 
away.  That  Is,  as  the  trail Ing-edge  velocity  changes  from  small  and  out¬ 
board  to  small  and  inboard  the  wake  position  a  finite  distance  downstream 
does  not  "flip  flop,"  but  changes  only  slightly.  The  question  is  how  can 
this  be  resolved  with  results  of  reference  11. 
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The  only  explanation  appears  to  be  that  as  the  trai ling-edge  velocity 
component  becomes  small  —  either  inboard  or  outboard  —  the  wake  approaches 
the  trailing  edge  bisector  at  small  finite  distances  from  the  trailing 
edge.  That  is,  the  curvature  of  the  wake  at  the  trailing  edge  becomes  large 
as  the  velocity  becomes  small  and  in  a  very  small  distance  the  wake  "curves 
around"  and  approaches  the  trailing-edge  bisector.  The  situation  is  sketched 
In  figure  14.  The  wake  approaches  the  trailing-edge  bisector  more  and  more 
rapidly  as  the  velocity  component  along  the  trailing  edge  approaches  zero. 

The  above  argument  requires  only  continuity  and  symmetry. 

Thus,  if  the  Kutta  condition  in  the  form  of  property  (a)  is  applied,  a 
finite  distance  from  the  trailing  edge  (as  it  must  be  in  the  present  method 

and  those  of  references  4,  5,  6  and  7)  and  if  the  sweep  angle  of  the 

trailing  edge  is  such  that  the  component  of  velocity  along  the  trailing  edge 
is  not  large,  then  the  point  must  lie  along  the  trailing-edge  bisector.  For 
example,  the  method  of  reference  6  applies  property  (a)  a  distance  of  3  percent 
of  local  airfoil  chord  downstream  from  the  trailing  edge  and  states  that  the 

point  should  lie  along  the  bisector  of  the  trailing  edge  rather  than  the 

tangent  to  the  upper  surface  as  required  by  reference  11.  On  the  other  hand. 


Figure  14.  Behavior  of  the  vortex  wake  near  the  trailing  edge  for 
small  values  of  the  trailing  edge  velocity  component. 
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It  is  clear  that  for  large  positive  sweep  angles,  the  component  of  velocity 
along  the  trailing  edge  becomes  the  same  order  of  magnitude  as  freestream 
velocity.  Presumably,  the  Kutta  condition  should  then  be  applied  along  the 
tangent  to  the  upper  surface.  At  what  value  of  trailing-edge  velocity  it 
becomes  necessary  to  change  from  one  scheme  to  the  other  is  not  known,  but 
it  certainly  must  be  dependent  on  the  distance  of  the  point  of  application 
from  the  trailing  edge.  Numerical  experiments  presented  in  reference  6  and 
similar  experiments  performed  by  the  present  author  show  that  the  calculated 
results  are  rather  sensitive  to  the  direction  of  the  point  of  application  of 
the  Kutta  condition.  For  a  typical  trailing-edge  angle  the  calculated  lift 
coefficient  obtained  from  an  application  point  on  the  trailing-edge  bisector 
can  easily  differ  by  20  percent  from  the  lift  coefficient  obtained  from  an 
application  point  on  the  upper-surface  tangent. 

Even  when  the  direction  from  the  trailing  edge  of  the  point  of  application 
of  the  Kutta  condition  is  not  a  problem,  calculations  using  property  (a) 

(wake  tangency)  are  affected  by  the  distance  of  the  point  of  application  from 
the  trailing  edge.  This  Is  to  be  expected.  What  is  not  expected,  however,  is 
that  If  property  (b)  (pressure  equality)  is  used  in  the  manner  described  above, 
the  calculated  results  are  not  sensitive  to  distance  from  the  trailing  edge. 

A  study  was  made  In  two-dimensional  flow,  where  the  streamline  is  known  to 
leave  the  trailing  edge  along  its  bisector.  The  airfoil  selected  was  a  sym¬ 
metric  10-percent  thick  RAE  101  airfoil  at  10  degrees  angle  of  attack.  Bound 
vortlcity  was  provided  by  a  constant-strength  sheet  of  vorticity  coincident 
with  the  airfoil  surface  as  described  In  section  6.3.  Cases  were  run  with 
27,  53,  and  103  surface  elements.  The  results  were  also  extrapolated  to 
infinite  element  number.  Calculated  lift  coefficients  are  shown  in  figure  15. 
Since  property  (b)  (pressure  equality)  is  applied  at  the  control  points  of 
the  two  elements  adjacent  to  the  trailing  edge,  there  is  just  one  lift 
coefficient  for  each  element  number.  These  are  plotted  at  the  chordwise 
distance  of  the  nearest  control  point  from  the  trailing  edge,  which  ranges 
from  1.75  percent  chord  for  the  27  element  case  to  0.25  percent  chord  for 
the  103  element  case.  Remarkably,  the  calculated  lift  coefficients  are  almost 
constant  at  a  value  of  about  0.944,  and  the  extrapolation  to  the  trailing 
edge  itself  (Infinite  element  number)  yields  a  value  of  0.943.  For  each 
element  number  property  (a)  (wake  tangency)  was  applied  along  the  trailing- 
edge  bisector  at  distances  from  the  trailing  edge  ranging  from  0.25  percent 
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ira  the  trailing  edge  of  the  point  of  application  of  the  Kut 
percent  thick  symmetric  section  at  10  degrees  angle  of  att 


chord  to  4.0  percent  chord.  The  calculated  lift  varies  significantly  with 
both  element  number  and  distance  of  the  application  point  from  the  trailing 
edge.  It  appears  that  results  are  more  sensitive  to  distance  from  the 
trailing  edge  than  to  element  number.  If  results  are  extrapolated  both  to 
infinite  element  number  and  to  zero  distance  from  the  trailing  edge,  the  lift 
coefficient  is  given  as  0.942.  This  agrees  with  the  value  obtained  by 
extrapolating  property  (b)  and  with  the  value  of  0.9423  obtained  by  a  high- 
accuracy  conformal  mapping  solution.  However,  for  the  27  element  case 
(a  reasonable  number  in  three  dimensions)  the  extrapolated  lift  coefficient 
for  zero  distance  is  0.955,  which  is  reasonably  close  to  the  correct  value, 
hut  use  of  a  point  of  application  at  3  percent  chord,  as  called  for  in 
reference  6,  gives  a  lift  coefficient  of  1.005,  which  is  considerably  in  error. 
Even  for  the  extrapolation  to  infinite  element  number,  a  point  of  application 
at  3-percent  chord  gives  a  lift  coefficient  of  0.989.  Thus,  it  appears  that 
use  of  a  pressure-equality  Kutta  condition  applied  on  the  body  (property  (b)) 
is  more  accurate  and  less  sensitive  than  the  flow-tangency  Kutta  condition 
applied  in  the  wake  (property  (a)),  which  is  used  in  references  4,  5,  6  and  7 
even  if  the  direction  by  which  the  wake  leaves  the  trailing  edge  is  not  a 
problem. 


6.6  Symmetry  Planes 

To  conserve  computing  time  and  reduce  the  required  input,  the  method  is 
equipped  to  take  advantage  of  any  planes  of  symmetry  the  flow  may  possess. 
Either  one  or  two  symmetry  planes  may  be  accounted  for.  The  xz-plane  is 
denoted  the  first  symmetry  plane.  If  there  is  one  plane  of  symmetry,1t  must 
be  the  xz-plane, and  the  points  iefining  the  nonredundant  portion  of  the  body 
must  be  input  accordingly.  The  xy-plane  is  denoted  the  second  symmetry  plane. 
If  there  are  two  symmetry  planes,  they  must  be  the  xz-  and  xy-planes,  and  the 
input  points  must  reflect  this.  Each  symmetry  plane  is  designated  either 
"plus"  or  "minus'.1  A  plus  symmetry  plane  has  zero  normal  velocity  at  all 
points  of  the  plane,  l.e.,  it  behaves  as  a  solid  wall.  A  minus  symmetry 
plane  has  zero  velocity  potential  at  all  points  of  the  plane.  The  usual 
application  in  aerodynamics  consists  solely  of  plus  symmetry  planes.  An 
example  of  a  negative  symmetry  plane  is  a  free  surface  for  the  condition  of 
infinite  Froude  number.  Thus,  a  hydrofoil  traveling  near  the  water  surface 
has  two  symmetry  planes  -  one  plus  and  one  minus. 
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Symmetry  Is  accounted  for  in  the  part  of  the  calculation  devoted  to 
the  velocity  induced  by  the  quadrilateral  surface  elements.  Recall  that  an 
element  may  have  on  It  either  a  source  or  a  dipole  distribution  or  both. 
Velocities  are  computed  at  all  control  points  due  to  the  source  and/or  dipole 
distribution  on  a  basic  element  defined  by  input  points.  Then  this  element 
Is  reflected  successively  in  the  symmetry  planes,  induced  source  and/or  dipole 
velocities  at  the  control  points  are  computed  for  the  reflected  elements, 
and  the  Induced  velocities  for  the  reflected  elements  are  added  to  the  cor¬ 
responding  quantities  for  the  basic  element.  Reflection  In  a  plus  symmetry 
plane  requires  a  source  distribution  of  the  same  sign  as  the  original  but 
a  dipole  distribution  of  opposite  sign.  (All  magnitudes  of  course  are  equal 
to  the  original.)  A  minus  symmetry  plane  yields  the  opposite  situation,  i.e., 
source  changes  sign,  dipole  does  not. 

In  symmetry  cases  it  is  assumed  that  the  y-direr.tion  is  essentially 
I  "spanwise"  on  the  wing,  so  that  the  first  symmetry  plane  is  the  midplane  of 

the  wing.  The  second  symmetry  plane  (if  any)  is  then  available,  for  example, 
as  a  ground  plane.  Figure  16  shows  a  section  of  an  element  and  its  bordering 
,  N-lines,  together  with  their  reflections.  The  N-lines  are  labeled  “first" 


Figure  16.  Reflection  of  an  element  and  its  associated  N-lines  in 
symmetry  planes. 
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and  "second'J  and  In  the  case  shown  the  "first"  N-llne  is  Inboard  from  the 
"second"  one  with  respect  to  the  span  direction  on  the  basic  element.  It  can 
be  seen  that  reflection  In  the  y-dlrectlon  reverses  this  relationship  while 
reflection  In  the  z-directlon  does  not.  This  condition  affects  the  assembly 
of  the  dipole  velocities,  and  thus  the  Input  points  should  be  compatible  with 
the  above  assumption. 


6.7  Multiple  Angles  of  Attack 

The  method  can  calculate  flow  about  a  lifting  configuration  for  several 
angles  of  attack  of  the  freestream  in  essentially  the  same  computing  time 
as  that  required  for  a  single  angle  of  attack.  In  the  latter  case,  sets  of 
source  density  are  obtained  for  L  +  1  onset  flows  —  1  uniform  stream  at 
angle  of  attack,  and  L  bound  vortlclty  onset  flows.  Here  L  is  the  number 
of  lifting  strips  and  Is  generally  in  the  range  10  to  30.  The  Kutta  condition 
then  yields  L  combination  constants  for  these  vorticity  flows.  It  is  also 
possible  to  Input  several  angles  of  attack,  say  F,  and  obtain  F  +  L  basic 
source  distributions  for  the  F  uniform  flows  and  L  vortlclty  flows. 

Then  the  Kutta  condition  is  applied  to  each  of  the  F  uniform  stream  solu¬ 
tions  separately  to  obtain  a  complete  set  of  L  combination  constants  for  the 
vorticity  flows.  Using  these  constants,  a  "combined"  source  density  distri¬ 
bution  is  obtained  for  each  angle  of  attack  in  the  manner  described  in  section 
6.5.  The  output  consists  of  a  complete  set  of  surface  pressures  and  off-body 
velocities  for  each  angle  of  attack.  For  comparison  purposes  nonlifting 
solutions  may  also  be  obtained  by  computing  strictly  from  the  source  densities 
obtained  from  the  uniform  onset  flows. 

When  computed  In  the  above  manner,  the  solutions  for  all  angles  of  attack 
have  the  same  position  for  the  trailing  vortex  wake.  This  Is,  of  course,  an 
approximation,  because  the  true  position  of  the  wake  changes  with  angle  of 
attack.  However,  as  will  be  discussed  In  section  8.5,  the  calculated  surface 
pressures  are  very  insensitive  to  angle  of  inclination  of  the  wake  with  the 
range  of  practical  interest.  Thus,  solutions  obtained  by  the  multiple  angle- 
of-attack  option  are  essentially  as  accurate  as  can  be  expected  from  potential 
flow. 
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6.8  Some  Special  Situations 


The  basic  theoretical  difficulties  with  the  potential -flow  model  for 
lift  (section  5.3)  have  their  effect  on  the  method  of  solution  by  necessitating 
special  handling  of  certain  frequently  occurring  situations.  The  special 
features  that  have  been  built  into  the  r. ethod  to  handle  these  situations  are 
discussed  In  this  and  In  the  following  section.  Other  special  situations  may 
be  discovered  In  the  future. 

Two  special  situations  exist  where  elements  must  be  placed  inside  the 
body  surface.  No  normal -velocity  boundary  condition  can  be  applied  at  such 
elements  and  no  source  density  should  be  applied  to  them.  Thus,  these  elements 
do  not  count  as  far  as  the  matrix  of  Induced  velocities  Is  concerned.  How¬ 
ever,  they  do  have  dipole  distributions  and  these  must  be  accounted  for. 

The  first  situation  occurs  when  a  nonlifting  portion  of  the  body  inter¬ 
sects  a  lifting  portion  at  a  finite  angle  (often  nearly  normal)  without 
breaking  the  continuity  of  the  trailing  edge.  An  example  is  provided  by 
the  wing-pylon  intersection  shown  in  figure  17.  A  certain  portion  of  the 


Figure  17.  Handling  of  a  wing-pylon  intersection. 
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lifting  body  surface  is  "Inside"  the  pylon.  Howeeer,  tfu-  d i ,  ,-v  i *<  distnouti 
should  be  continuous  through  this  reqion  to  ;• . v o  1  < i  numerical  difficulties. 
Thus,  as  far  as  dipole  calculations  are  cono-mmi,  *he  "inside"  e Insert ^  c,r 
normal  members  of  the  liftinq  strips  to  which  they  belong.  But  they  are- 
ignored  as  far  as  source  calculations  or  bc-.niar/  conditions  arc  c orcet  ned . 
Such  elements  are  designated  "ignored  elements They  usually  conga  i;,e  or1 
part  of  a  lifting  strip. 

The  second  situation  occurs  when  a  lifting  portion  of  the  body  intense 
a  nonlifting  portion  at  a  finite  angle  (often  ne^-  i  y  n  .  T ►. i»  inu.no  tar 
case  of  this  Is  the  wing-fuselage  intersection,  as  ii  •u'-4ra<ed  in  tlou-r  U 
As  Is  well-known,  the  locaT'section  lift  coefficient"  c-n  the  wing  does  not 


Figure  18.  Handling  of  a  wing-fu-.eiogr  sr'f-sec 
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fall  to  zero  at  the  fuselage  intersection.  Thus  the  dipole  strength  on  the 
N-Hne  lying  along  the  Intersection  is  not  zero.  However,  the  lifting  section 
cannot  simply  be  terminated,  because  that  would  result  in  a  concentrated 
vortex  filament  right  on  the  surface.  Accordingly,  an  additional  or  "extra" 
lifting  strip  Is  added  to  the  lifting  section  (see  figure  18).  It  is  either 
the  first  or  the  last  strip  of  the  lifting  section.  The  extra  strip  lies 
inside  the  nonlifting  body  and  is  a  complete  lifting  strip  including  wake. 

No  source  densities  or  normal-velocity  boundary  conditions  are  applied  to 
the  elements  of  the  extra  strip.  The  dipole  strength  is  taken  constant  in 
the  "spanwlse"  direction  across  the  extra  strip.  The  value  of  the  dipole 
strength  on  the  extra  strip  is  such  that  the  dipole  strength  is  continuous 
across  the  N-line  lying  along  the  wing-fuselage  intersection.  The  interior 
edge  of  the  extra  strip  has  nonzero  dipole  strength  and  may  lead  to  a  concen¬ 
trated  vortex  In  the  streamwiso  direction.  For  example,  shown  in  figure  18, 
the  vortex  may  lie  along  the  fuselage  centerline  and  its  downstream  extension. 
If  the  lifting  configuration  has  a  right-and-left  symmetry,  e.g.,  a  fuselage 
with  both  wings,  and  i.  ♦he  flow  is  also  symmetric,  e.g.,  zero  yaw,  the 
extra  strips  for  the  right  and  left  sides  have  the  same  strengths  on  their 
Interior  edges.  Tnur,  in  this  case  there  is  no  discontinuity  of  dipole 
strength  and  no  concentrated  vortex.  If,  however,  the  lift  Is  not  symmetric, 
there  will  be  a  concentrated  vortex.  This  Is  unavoidable  Decause  ;t  is 
physically  real.  An  example  is  the  hub  vortex  or  a  propeller.  This  also 
occurs  at  a  tip  tank,  which  is  essentially  a  small  fuselage  with  only  one 
wing.  However,  the  case  shown  in  section  1C,1  exhibits  no  numerical  difficulty. 

6,9  Summary  of  the  Logic  of  the  Calculation 

The  overall  logic  of  the  method  is  rather  similar  to  that  of  the  method 
for  nonlifting  potential  flow  described  in  section  6.1.  There  are,  of  course, 
certain  additions.  The  order  of  the  various  parts  of  the  calculation  and 
their  functions  are  outline''!  below. 

The  geometry  of  the  three-dimensional  configuration  is  input  to  the 
program  In  the  form  of  coordinates  of  a  set  of  points.  The  points  are  input 
along  N-lines,  which  are  essentially  coordinate  curves  in  the  body  surface 
(figure  G).  The  configuration  is  divided  into  lifting  and  nonlifting  portions 
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as  discussed  in  section  6.2.  Each  of  these  may  be  further  divided  into 
sections  -  lifting  and  nonlifting.  The  nonlifting  sections  are  input  first. 

The  N-lines  of  the  lifting  section  define  both  the  body  and  the  trailing 
vortex  wake.  Coordinates  of  points  off  the  body  in  the  flow  field  where 
flow  calculations  are  desired  are  Input  after  the  points  defining  the  body. 

If  the  Kutta  condition  is  applied  by  a  condition  of  flow  tangency  downstream 
of  the  trailing  edge  in  the  wake,  coordinates  of  these  points  and  the  cor¬ 
responding  normal  vectors  are  Input  between  the  on-body  and  the  off-body 
points.  The  remaining  input  consists  of  control  flags  governing  the  logic 
of  the  calculation  plus  a  few  parameters,  such  as  angle  of  attack. 

Surface  elements  are  formed  from  input  points  in  the  manner  described 
in  section  7.2  for  lifting  sections  and  In  the  manner  of  reference  3  for 
nonlifting  sections.  The  "formation"  of  an  element  consists  of  the  calcula¬ 
tion  of  various  geometric  quantities  associated  with  that  element,  including 
coordinates  of  the  control  point  (centroid),  components  of  unit  vectors  along 
the  axes  of  a  coordinate  system  based  cn  the  element,  one  of  which  is  the 
unit  normal  vector  to  the  element,  and  moments  of  the  area  of  the  element. 
Elements  of  lifting  sections  are  logically  associated  into  lifting  strips 
of  elements,  which  consist  of  those  elements  formed  from  the  same  two  N-lines. 

Every  element  has  on  it  a  constant  source  density.  Lifting  elements  also 
have  a  dipole  distribution.  Formulas  have  been  derived  that  enable  velocities 
induced  by  the  elements  at  points  In  space  to  be  calculated  (section  7.0  for 
lifting  elements  and  reference  3  for  nonlifting  elements).  For  each  element 
the  velocities  induced  by  Its  constant  source  density  at  all  control  points 
and  off-body  points  are  computed  and  saved  in  low-speed  storage  (tape  or  disk). 
If  there  are  symmetry  planes,  velocities  induced  by  the  reflections  of  an 
element  are  added  to  the  velocities  due  to  the  element  itself.  This  is  the 
vector  matrix  of  induced  velocities.  For  each  element  of  a  lifting  section 
velocities  Induced  by  Its  dipole  distribution  at  all  control  points  and 
off-body  points  are  computed.  These,  however,  are  not  saved  Individually. 
Instead,  dipole  velocities  for  all  elements  of  a  lifting  strip,  including 
wake  elements,  are  added  to  obtain  velocities  due  to  the  entire  strip.  Thus, 
if  there  are  N  source  elements,  at  whose  control  points  the  normal-velocity 
boundary  condition  Is  to  be  applied,  0  off-body  points,  and  L  lifting 
strips,  there  is  a  (N  +  0)  x  N  matrix  of  induced  source  velocities  and  a 
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(N  +  0)  x  L  matrix  of  induced  dipole  velocities,  for  the  "step  function" 
bound  vorticity  option,  the  dipole  (vorticity)  velocities  induced  by  a  lifting 
strip  are  those  due  to  a  spanwise  constant  dipole  distribution  and  they  can 
be  computed  in  a  straightforward  manner.  For  the  "piecewise  linear"  bound 
vorticity  option,  two  sets  of  Induced  dipole  velocities  are  computed  for 
each  lifting  strip:  one  due  to  a  spanwise  constant  dipole  distribution  and 
one  due  to  a  linear  distribution  with  unit  rate  of  change  in  the  spanwise 
direction  and  zero  value  at  the  "midspan"  of  the  strip.  These  are  then 
combined  using  the  mechanism  of  the  parabolic  fit  and  the  conditions  at 
the  ends  of  the  lifting  sections  to  obtain  L  sets  of  induced  dipole 
velocities,  each  of  which  is  proportional  to  the  midspan  value  of  bound 
vorticity  on  one  lifting  strip.  The  calculations  outlined  in  this  paragraph 
comprise  one  of  the  two  time-consuming  parts  of  the  method. 

The  first  N  rows  of  the  induced  source  velocity  matrix  are  the 
velocities  at  the  control  points.  Components  of  these  velocities  along  the 
local  normal  direction  are  computed  to  yield  an  N  x  N  scalar  matrix  of 
Induced  normal  velocities.  This  is  the  coefficient  matrix  of  the  linear 
equations  for  the  source  density.  The  right  sides  of  the  linear  equations 
consist  of  components  along  the  local  normal  direction  of:  F  uniform  onset 
flows  at  various  angles  of  attack  and  the  first  N  rows  of  the  induced 
dipole  velocity  m^crix.  The  linear  equations  are  solved  by  direct  elimination 
to  yield  (F  +  L)  sets  of  source  densities  on  the  N  source  elements. 

The  matrix  solution  is  the  second  time-consuming  part  of  the  method. 

Flow  velocities  are  computed  for  all  (F  +  L)  sets  of  source  density 
at  the  points  used  to  establish  the  Kutta  condition.  These  are  the  2L 
control  points  adjacent  to  the  trailing  edge  on  all  lifting  strips  if  the 
condition  of  equal  upper  and  lower  surface  pressure  is  used.  If  the  condition 
of  flow  tangency  in  the  wake  downstream  of  the  trailing  edge  is  used,  the 
points  are  L  particular  off-body  points  Input  to  the  program.  The  Kutta 
condition  is  formulated  as  L  equations  for  the  L  values  of  bound  vorticity 
on  the  lifting  strips  using  each  of  the  F  uni  form- stream  flows  in  turn 
with  the  L  vorticity  flows.  The  result  is  F  sets  of  L  values  of  bound 
vorticity.  For  each  uniform  onset  flow  a  "combined"  set  of  source  densities 
Is  obtained  as  a  linear  combination  of  the  basic  L  sets  of  source  densities 
corresponding  to  the  vorticity  flows  and  the  set  of  source  densities  for 
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the  uniform  flow  itself.  The  combination  constants  for  the  vorticity  flows 
are  values  of  bound  vorticity  obtained  from  the  Kutta  condition.  There  are 
F  sets  of  N  values  of  the  combined  source  densities.  Similarly,  the  same 
values  of  bound  vorticity  are  used  as  combination  constants  to  obtain  a 
"combined"  onset  flow  at  all  N  +  0  points  where  velocities  are  to  be 
computed.  There  are,  of  course,  F  such  "combined"  onset  flows. 

A  complete  flow  solution  is  computed  using  each  set  of  combined  source 
densities  and  its  corresponding  combined  onset  flow.  Such  a  solution  consists 
of:  flow  velocities  and  pressures  at  all  control  points,  flow  velocities  at 
all  off-body  points,  the  bound  vorticity  values  used  to  satisfy  the  Kutta 
condition,  and  integrated  forces  and  moments  on  each  lifting  strip,  on  each 
lifting  and  nonlifting  section,  and  on  the  entire  configuration.  An  option 
also  exists  for  computing  a  nonlifting  solution  at  each  angle  of  attack  by 
setting  all  values  of  bound  vorticity  equal  to  zero. 


'  '  '  1 . ' . '  ; . . 1  . . m -li , a ; -1 


7.0  DETAILS  OF  THE  METHOD  OF  SOLUTION 


7.1  Order  of  the  Input  Points 

As  mentioned  previously,  the  points  defining  the  body  surface  are  input 
N-line  by  N-iine,  and  the  points  on  a  given  N-line  are  input  consecutively. 

The  order  of  the  input  determines  the  direction  of  the  outer  normal  vectors 
to  the  elements,  i.e.,  determines  whether  the  case  in  question  is  an  interior 
or  an  exterior  flow.  The  rule  for  insuring  that  normal  vectors  point  into 
the  field  of  flow  rather  than  into  the  interior  of  the  body  is  the  same  as  in 
reference  3:  If  an  observer  in  the  field  of  flow  looking  towards  the  body 
surface  sees  N-lines  input  from  left  to  right,  he  should  also  see  individual 
points  on  an  N-line  input  from  bottom  to  top.  An  example  of  correct  input  for 
the  right  wing  of  an  airplane  is  as  follows:  The  N-lines  are  input  from  tip 
to  root.  On  each  N-line  the  points  are  input  beginning  at  the  trailing  edge, 
traversing  the  lower  surface  to  the  leading  edge,  returning  to  the  trailing 
edge  along  the  upper  surface,  and  continuing  into  the  wake.  The  alternate 
way  of  inputting  a  right  wing  is  to  input  the  N-lines  from  root  to  tip  and 
on  each  N-line  to  input  upper  surface  points  first  followed  by  the  lower 
surface  points  and  the  wake  points.  Both  of  these  input  schemes  produce 
identical  surface  elements.  However,  they  lead  to  somewhat  different 
implied  surface  dipole  distributions.  This  matter  is  discussed  in  section  7.3. 
The  conclusion  is  that  the  first  of  the  two  input  schemes  above  is  to  be 
preferred.  In  any  case,  the  logic  of  the  program  for  determining  which 
elements  are  on  the  surface  and  which  are  on  the  wake  requires  that  the  first 
point  on  an  N-line  of  a  lifting  section  be  at  the  trailing  edge. 

7.2  Formation  of  the  Elements  from  Input  Points 

This  section  outlines  the  way  that  the  elements  are  actually  formed 
from  the  input  points.  There  are  two  principal  differences  between  the 
formation  of  lifting  elements  and  that  of  nonlifting  elements.  The  first 
Is  the  manner  of  adjusting  the  input  points  to  make  a  plane  quadrilateral. 

The  second  is  the  calculation  of  area  moments  up  to  fourth  order.  The 
procedure  for  forming  nonlifting  elements  is  given  in  reference  3  and  will 
not  be  repeated  In  this  section,  which  is  concerned  solely  with  lifting 
elements . 
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Let  the  coordinates  of  the  input  points  that  are  used  to  form  an 
1  i  1 

element  be  denoted  x^,  y^,  z^,  k  =  1,  2,  3,  4.  These  coordinates  are  with 
respect  to  the  reference  coordinate  system,  the  system  in  which  the  physical 
lifting  configuration  is  defined.  It  simplifies  the  equations  to  use  vector 
notation,  so  define 

V  xk  ^  +  y!  ^  +  z!  ^  (7.2.1) 

where  T,  J,  f,  are  unit  vectors  along  the  axes  of  the  reference  coordinate 
system.  Recall  that  an  element  is  formed  from  points  on  two  consecutive 
N-lines.  The  input  points  k  =  1  and  2  are  on  one  N-line,  the  "first" 

N-line,  and  the  Input  points  k  =  3  and  4  are  on  the  "second"  N-line.  In 
what  follows,  subscripts  F  and  S  are  used  to  denote  quantities  associated 
with  the  first  and  second  N-lines,  respectively.  The  numbering  k  =  1,  2,  3,  4 
is  cyclic  around  the  element  to  be  consistent  with  reference  3.  The  adjustment 
of  the  input  points,  which  is  shown  in  figure  IS,  is  as  follows. 

First  form  the  N-line  vectors 


figure  19.  Adjustment  cf  the  input  points  to  form  a  plane  trapezoidal 
element. 
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The  two  parallel  sides  of  the  trapezoid  are  taken  as  parallel  to  the  weighted 
average  of  these  two  vectors.  In  the  coordinate  system  of  the  element  this 
is  also  the  direction  of  the  x-axis.  The  unit  vector  parallel  to  the  two 
parallel  sides  of  the  trapezoid  is  denoted  to  show  It  is  also  the  unit 

vector  along  the  x  or  £  axis  of  the  element  coordinate  system.  It  is 
computed  from 


'E  = 


p  +  f 

kf  ks 


(7.2.3) 


where  |v|  means  absolute  magnitude  of  the  vector  v,  i.e.,  the  square  root 
of  the  sum  of  the  squares  of  its  components.  The  parallel  sides  have  the 
direction  of  T^.  The  calculation  insures  that  each  parallel  side  has  the 
same  midpoint  and  the  same  length  as  the  segment  of  N-line  from  which  it  is 
formed.  In  fact,  once  the  elements  are  formed  the  original  N-line 
segments  are  replaced  by  these  parallel  sides.  The  side  lengths  are 


dF=  ?f 


V  Fs 


(7.2.4) 


The  midpoints  in  vector  form  are 


Xp  -  ■jr  (  x  .j  +  X  j ) , 


x$  =  (*3  +  *4) 


(7.2.5) 


The  endpoints  of  the  two  parallel  sides,  which  are  thus  the  corner  points  of 
the  trapezoidal  element  are,  in  vector  form, 


x  -I  —  X  r~  —  yr  d  . 

1  r  2  F  E 


x3  =  XS  +  1  'Ve  ’ 


+  i  dcTr 


X2  =  XF  T  7  UF1E 


x4  ”  XS  7  ds’E 


(7.2.6) 


The  normal  vector  to  the  plane  of  the  elemi  nt  is 


N  *  (x4  -  x2)  x  (x3  -  x] ) 


(7.2.7) 


The  unit  normal  vector  is 


n  =  JL 

INI 


(7.2.8) 
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This  Is  also  the  unit  vector  along  the  z-axis  of  the  element  coordinate  system. 
The  unit  vector  along  the  y  or  n  axis  of  the  element  coordinate  system  Is 

TE  =  nxTE  (7.2.9) 

In  comp  xnent  form  the  three  unit  vectors  are 

iE  =  +  ai2^"  +  a13^ 

Te  =  a2l^  +  a22^  +  a23^  (7.2.10) 

n  =  kE  '  u31^  +  a32^"  +  a33^ 

The  3  x  3  array  of  a's  is  the  transformation  matrix  that  is  used  to  trans¬ 
form  coordinates  of  oolnts  and  components  of  vectors  between  the  reference 
and  element  coordinate  systems  In  the  manner  described  in  reference  3. 


Temporarily  the  origin  of  the  element  coc  dlnate  system  is  taken  as  the 
point  whose  coordinates  are  the  averages  of  those  of  the  Input  points.  (The 
same  averages  are  obtained  using  the  corner  points.)  In  vector  notation,  the 
average  point  Is 

With  this  origin,  the  element  coordinates  of  the  corner  points  are 


«k  '  aH  (xk  -  xav>  *  a12^yk  -  yj  +  a13(zk  -  zav> 

nk  ’  a21(xk  ~  xa*'  +  a22^yk  -  yav*  *  a23^zk  _  zav* 
k  =  1,2, 3, 4 


(7.2.12) 


where  In  accordance  with  vector  notation,  x^,  z^  are  the  coordinates  of 
from  (7.2,6).  It  will  turn  out  that 


*  * 

n-|  ”  n2  and 


(7.2.13) 
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The  width  of  the  element  is 


★  *  ir 

w  =  n,  -  n3  •  2n1 


The  slopes  of  the  nonvertical  sides  of  the  element  (figure  20)  are 


e2“*3 


_  *1  -*4 


'"32  "  w  41  '  ~ w 

with  respect  to  the  m  axis.  The  coordinates  of  the  centroid  are 


w 


m32  m41 


'o  6  *  *  *  * 

c3  +  €2  -  h  -  e4 


5o  = 


m32  +  m41 


~2 - no 

The  reference  coordinates  of  the  centroid  are 


xo  =  xav  +  all£o  +  a21 no 


yo  =  yav  +  a12^o  +  a22no 


zo  =  Zav  +  a13f>o  +  a23no 


Figure  20.  A  plane  trapezoidal  element. 


The  centroid  is  now  taken  as  the  origin  of  the  element  coordinate  system  and 
replaces  the  average  point  in  all  subsequent  calculations.  With  respect  to 
the  centroid  as  origin,  the  element  coordinates  of  the  corner  points  are 


^k  "  ^k 


nk  =  \  ~  no 


(7.2.18) 


where 


n2  =  nl 


n4  =  n3 


(7.2.19) 


These  are  the  corner  point  coordinates  used  In  all  subsequent  calculations. 

Several  other  geometric  quantities  are  needed  in  future  calculations. 
These  are  now  computed.  The  Intercepts  where  the  sides  intersect  the  x  or 
i  axis  (figure  20)  are 


^3n2  ^2n3 

b32  =  - w - 


_  *>4nl  ~  *hn4 


(7.2.20) 


The  maximum  diagonal  of  the  element  is 


t  =  Max 


Ug  ~  S4)  +  (ig  —  n^)‘ 
^U3  -  S])Z  +  (n3  - 


(7.2.21) 


The  lengths  of  the  sides  are 


d12  =  dF 


d32  =  +  m32 


d34  =  dS 


d41  =  +  m4i 


(7.2.22) 


Also  needed  are  the  total  arc  lengths  along  the  N-lines  from  the  trailing  edge 
up  to  the  element  in  question.  These  are 


•F  =  Zdr*  LS  =  Yj' 


(7.2.23) 


where  the  sums  are  over  previous  elements  of  the  lifting  strip. 
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Finally  the  normalized  moments  of  the  area  of  the  element  are  required. 
These  are  defined  by 


U  ■  ^  <7-2-24> 

E 

where  the  region  of  integration  is  the  area  of  the  element.  For  example, 

2  4  4  4 

t  Iqq  is  the  area  of  the  element,  t  t  1^,  t  Iq2  are  the  moments  of 
Inertia  or  second  moments.  The  first  moments  and  *01  are  zero  because 

the  centroid  is  used  as  origin  of  coordinates.  The  order  of  a  moment  is  the 
sum  of  its  subscripts  n  +  m.  There  are  three  second-order  moments,  four 
third-order,  and  five  four' h-order.  The  present  method  uses  up  through  fourth 
order.  The  moments  are  calculated  by  a  straightfc  ward  but  rather  lengthy  set 
of  formulas  which  are  given  below. 

First,  normalize  the  corner  point  coordinates  by  the  maximum  diagonal, 

=  Sk/t,  \  =  k  =  1,  2,  3,  4  (7.2.25) 

Now  the  normalized  moment  may  be  defined  in  terms  of  certain  auxiliary 
functions 


,  _  t(32)  .  T (41 )  .  1  r-m+l,:n+l  :n+lx  ,  -m+l,;n+l  ;n+lx-, 

!nm  "  ”*nm  +  V,  *  (iTT W’^TT  [nl  U2  '  ?1  }  +  r,3  (f4  ~  ?3  )] 

(7.2.26) 


The  auxiliary  function  I 
If  |m32|  >  1: 


(32) 


nm 


is  as  follows1 


.(32)  _  1  rZn-H  ’m+1-,2 

nm  "  (m  +  1  )(n  +  1 )  *■ '  n  ■'3 

1  1  r’n+2  *m-i2 

(n  +  1  j (n  +  2)  mjg  ?  n  j 


*  TFTrmrrrnTnTT  x  u"+3  .■'h? 

1 


m32  "  "3 


(7.2.27) 


fil(m  1  )  1  •  ^j./i  o  fy 

-  Tn  - TTTn  +  : 2)ln  +  3T(n”+  4)  Un+4  nm_2]2 
f  m(m  *“  1 ) (m  2)  1  •,*,  o  o 

*  X  TRW V  TTCn  +  3 Hn  *  4)(n  +  5)  T  U"*5  n”"3]? 

32 

m(m  —  l)(m  —  2)(m  -  3)  1  .„lC  .„  A  o 

~  In  T  I:(n  T  2)(n_+ mf^7)Tn~5TTfr+“5T  ^  U"  6  n"1'4]^ 
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If  Im 


"32  - 


<  1: 


,(32)  _  1  _  r!n  'm+2n2 

!nm  '  Tm_+  \)(m  +  T)  m32  ^  n  J3 

n  _2  rln-1  *m+3T2 

~  (m  +  1‘) (in  “  ?)(m  +  3)  m32  U  n  ]3 

,  n(n  —  1)  _3  r>-2  >+4,2 

Tin  +  1 ) ( m  +  7J (m  +  31  (m  V  47  32  Lc  n  ]3 


(7.2.28) 


n(n  —  1 ) (n  —  2)  _4  rZn-3  >+5,2 

”  (m~Tr)"(“+  2)(iT+  3)Tm'  T  4)T'm  +  T)  m32  U  n  J3 

n(n  -  l)(n  —  2)(n  —  3)  _5  rIn-4  *m+6,2 

+  Tin  +  1  )(rii  +"  2')(m  +  ^{m  +  4Vm  +  5}{m  '+  €)'  32  U  n  J3 


where  the  bracketed  symbols  are  defined  by 

n*3]2  =  ^2  n2  —  ?3  n3  (7.2.29) 

(The  superscripts  in  the  above  equations  denote  simple  powers  of  the  quantities 
except  for  the  bracketed  double  superscript  (32),  which  denotes  the  side  of 
the  quadrilateral.)  It  is  clear  from  the  above  that  the  calculation  of  I^2^ 
requires  m  +  2  terms  of  (7.2.27)  or  n  +  1  terms  of  (7.2.28).  The  calcu¬ 
lation  is  simply  terminated  at  this  number  of  terms.  The  auxiliary  function 
(41 ) 

Inm  1S  stained  from  the  above  by  an  obvious  substitution  of  subscripts. 

All  the  above  geometric  quantities  associated  with  a  given  element  are 
saved  and  used  as  needed  to  calculate  velocities  induced  by  that  element. 

At  this  stage,  some  of  the  generated  quantities  are  output,  and  the  calcu¬ 
lation  may,  if  desired,  be  terminated.  The  purpose  of  this  option  is  to 
provide  an  opportunity  to  discover  errors  in  the  input  points  before  the 
execution  of  the  lengthy  calculations  of  the  main  part  of  the  program. 

7.3  Form  of  the  Surface  Dipole  Distribution 


7.3.1  General  Form.  Order  of  the  Input  Points 

The  surface  of  the  lifting  section  is  imagined  covered  with  a  dipole 
distribution  that  varies  in  the  following  manner.  The  dipole  strength  ^  is 
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fixed  as  zero  at  the  first  point  of  each  N-line  on  the  trailing  edge.  Along 
each  N-line  the  dipole  strength  is  proportional  to  distance  along  the  section 
curve.  This  curve  goes  completely  around  the  body  and  back  to  the  trailing 
edge,  at  which  point  p.  has  some  final  nonzero  value.  Behind  t.he  trailing 
edge  p  Is  constant  and  equal  to  Its  final  trailing  edge  value.  In  the  first 
input  example  of  secton  7.1,  a  right  wing  is  input  from  tip  to  root,  and 
points  on  an  N-line  traverse  the  lower  surface  to  the  leading  edge  and  return 
to  the  trailing  edge  along  the  upper  surface.  For  this  example  the  dipole 
variation  Is  as  shown  In  figure  21.  The  constant  of  proportionality  that 
expresses  the  variation  of  p.  along  each  N-line  Is  initially  unknown  and  its 
value  is  ultimately  determined  from  the  Kutta  condition.  Since  the  N-lines 
are  roughly  "chordwlse"  or  "streamwise"  on  the  lifting  surface,  this  constant 
Is  the  derivative  of  p  with  respect  to  distance  in  the  chordwise  direction. 
Thus,  according  to  the  result  of  Appendix  A,  the  proportionality  constant 
that  determines  the  growth  of  p  along  each  N-line  is  essentially  the 
"spanwise"  vorticity  strength  at  that  N-line,  which  is  the  bound  vorticity 
that  gives  the  lift. 

As  mentloneu  in  section  7.1,  points  along  an  N-line  of  a  lifting  section 
are  input  beginning  with  the  trailing  edge,  traversing  the  section  curve  of 
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Figure  21.  Variation  of  dipole  strength  a’ong  an  N-line. 
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the  wing  and  continuing  intc  the  wake.  The  orde»  may  be  sn-w 

so  that  either  the  lower  surface  is  Input  first,  as  i . lo  ;t rated  \n  ti.vn;  :-l, 
or  else  the  upper  surface  is  input  first,  in  which  ■  a s t-  rf-e  metier.  i . 

traversed  in  the  opposite  sense  to  that  shown  ir  figu'e  f '  (com  t*  <.  U.- w  m  ,' 
instead  of  clockwise).  Thus,  the  positive  direction  c .«  ace  ifnc  l>  a  the 
section  curve  is  opposite  in  the  two  input  scheres.  ’♦  the  oh-tc  >.*■  bound 
and  trailing  vorticity  are  to  be  Identical,  as  Mey  ■-•u-'it  ’>?  if  the  sa* v  m«,ly 
is  input  two  ways,  then  the  constants  of  proportional}  ‘>M t  :Ji}  *•'.*> 

strength  to  arc  length  along  the  N-line  for  each  of  t.v  t>.>  input  t chews 
must  turn  out  to  be  equal  in  magnitude  and  opposite  i -  : ign  To  V-Vstra?* 
the  two  cases,  let  the  constant  of  proportionality  t-e  :Vnr.»„oc  t-  .-’■■■••  an 
length  along  the  section  curve  be  denoted  s.  1  not.  ...  :  f:s  i:  th>  i 
order  is  that  of  figure  21,  and  p  =  -6s  if  tht  or  dr  'S  Toe 

dipole  strengths  along  the  N-line  for  the  two  c  cos  3re  i  i  Ivo-ruL*.-;  'v. 
figures  22a  and  22b.  The  dipole  strengths  in  th-.»  wakr  ?qjel  became  t‘<* 
reversal  of  sign  of  the  normal  vector  cancels  ot f  f«r  ••?rarc*<-«.  s>qr  reverse's 
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Figure  22.  Three  variations  of  dipole  stre-.gth  along  a  sect,--'  curve 

(N-line).  (a)  Clockwise  order  of  i-'p.,t  points  Hewer  su-  fate 
input  first),  (b)  Counterclock- ise  order  of  input  pr>i:  s 
(upper  surface  input  first).  )  C-jrs  tat*  ci  e  s  Irene  to 
on  body  obtained  by  subtracting  costs  (o'  ana  (t,. 
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of  the  dipole  strength  The  two  solutions  represented  by  figures  22a  and 
22b  may  be  subtracted  to  give  the  solution  illustrated  in  figure  22c,  which 
has  a  constant  dipoie  strength  B  •  s(tot.)  all  around  the  profile  curve  and 
zero  strength  in  the  wake.  Since  both  the  bound  vorticit.y  strength  B  and 
the  total  arc  length  around  a  section  curve  vary  with  "spanwise"  location, 
the  dipole  strength  of  figure  22c  varies  in  the  "spanwise"  direction  but  not 
in  the  "chordwise"  direction.  Thus,  according  to  Appendix  A,  the  related 
vorticlty  distribution  consists  of  closed  constant-strength  filaments  lying 
around  the  section  curves.  For  the  usual  case  of  a  wing  with  right  and  left 
symmetry  at  zero  yaw,  the  symmetry  insures  a  zero  total  strength  for  this 
vorticity  distribution.  Moreover,  the  flow  solution  of  figure  22c  has  no 
uniform  onset  flow,  which  was  canceled  in  the  subtraction  of  the  solutions 
corresponding  to  figures  22a  and  22b.  The  solution  corresponding  to  figure  22c 
is  continuous,  because  the  wake  singularity  is  zero,  and  It  satisfies  the 
classical  problem  defined  by  equations  (5.1.3),  (5.1.5),  and  (5.1.4)  with  zero 
right  side.  This  problem  has  a  unique  solution,  namely  the  trivial  solution. 
Thus,  the  solution  of  figure  22c  represents  zero  flow,  and  the  solutions  of 
figures  22a  and  22b  are  Identical  as  they  should  be. 

The  theoretical  considerations  of  the  previous  paragraph  are  strictly 
true  for  closed  bodies  in  the  limit  of  an  infinite  number  of  surface  elements. 
An  "open"  wing  tip  of  the  type  illustrated  in  figure  4c  is  excluded.  For 
practical  element  numbers,  numerical  experiments  must  be  performed.  Results 
of  such  an  experiment  are  presented  in  section  8.4  and  are  anticipated  here. 
When  a  wing  was  Input  in  the  two  ways  discussed  above,  the  resulting  bound 
vorticity  distributions  were  identical.  The  resulting  surface  pressure 
distributions  were  arly  identical  except  near  the  wing  tip  where  the  input 
scheme  illustrated  in  figures  21  and  22a  seemed  to  give  the  more  reasonable 
solution.  Accordingly,  it  was  concluded  that  points  on  N-llnes  of  lifting 
sections  should  be  input  with  the  lower  surface  first,  as  shown  in  figure  21. 
However,  two  wing-fuselages  input  with  the  upper  surface  first  have  very 
reasonable  surface  pressures.  The  preceeding  applies  to  positive  angle  of 
attack,  for  which  the  lower  surface  faces  the  onset  flow.  For  more  general 
flows  the  word  "lower"  in  the  above  is  replaced  by  "windward". 


7.3.2  Variation  Across  the  Span  of  a  Lifting  Stri 


The  variation  of  ^  between  the  two  N-llnes  used  to  form  a  lifting  strip 
Is  assumed  to  be  one  of  two  forms  that  correspond  to  the  "step  function"  and 
"piecewise  linear"  options  for  the  spanwise  bound  vorticity  variation,  as 
discussed  In  section  6.3.  For  the  "step  function"  option  the  proportionality 
constants  on  the  two  N-lines  bounding  the  strip  are  set  equal.  This  common 
value  Is  essentially  the  bound  vorticity  on  the  strip  and  is  determined 
directly  by  the  Kutta  condition.  In  general,  the  bound  vorticity  is  different 
on  adjacent  lifting  strips.  Thus,  there  are  really  two  values  of  "the"  pro¬ 
portionality  constant  on  an  N-line,  namely  those  corresponding  to  the  two 
lifting  strips  on  either  side  of  the  N-line.  The  dipole  distribution  is 
discontinuous  across  the  N-line,  which  Implies  a  discontinuity  of  bound  vor¬ 
ticity  and  a  concentrated  trailing  vortex  filament  along  the  N-line.  The 
"piecewise  linear"  option  essentially  assumes  a  linear  "spanwise"  variation 
across  a  lifting  strip  for  the  "chordwise"  proportionality  constant  of  the 
dipole  strength.  The  "spanwise"  derivative  is  determined  by  the  parabolic 
fit  discussed  In  section  7.11.  The  discontinuity  at  the  N-line  is  reduced 
to  a  higher  order  effect.  As  is  shown  below,  this  option  require'  an 
additional  dipole  term  in  the  wake. 


7.3.3  Variation  Over  a  Trapezoidal  Element 


Consider  now  a  typical  trapezoidal  lifting  element,  as  shown  in  figure  20. 
As  defined  in  section  7.2,  the  lifting  strip  to  which  the  element  belongs  Is 
bounded  by  two  N-llnes,  which  are  designated  the  "first"  N-line  and  the  "second" 
N-line  and  which  are  represented  by  subscripts  F  and  S,  respectively.  The 
constants  of  proportionality  for  th?  dipole  strength  along  the  N-llnes  are 
Dp  and  B<~,  respectively.  Thus,  if  s  denotes  arc  length  along  an  N-line: 


k  =  BpS  along  the  first  N-line 

h  =  B<.s  along  the  second  N-line 


(7.3.1) 


On  the  element  Itself,  the  parallel  side  at  n  =  is  a  segment  of  the 
first  N-line,  and  the  parallel  side  at  n  =  n3  Is  a  segment  of  the  second 
N-line  (figure  20).  The  dipole  strength  varies  linearly  along  each  side, 
namely 
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H  =  Ap 

+  BpS 

on 

n  =  Hi 

H  =  A$ 

+  B$C 

on 

n  =  n3 

(7.3.2) 


On  the  element  the  arc  length  along  the  N-lines  is  simply  the  coordinate  ? 
and  the  direction  of  increasing  arc  length  is  the  positive  s  direction.  On 
each  side  the  constant  A  is  the  value  of  n  for  c  =  0.  Thus, 


Ap  =  Bp  (tctal  arc  length  of  axis  from  trailing  edge) 


From  figure  20  and  equation  (7.2.23) 

hp  =  Lp  -  C1 

Similarly 


v/here 


hS  =  LS 


(7.3.3) 


(7.3.4) 


(7.3.5) 

(7.3.6) 


Now  the  dipole  distribution  K  on  the  element  is  assumed  in  the  form  of 
a  general  two-variable  second  degree  polynomial.  When  conditions  (7.3.2)  are 
applied,  it  turns  out  that  n  must  have  the  form 


Bc  —  Bc  Ac  —  Ac  Bcn-i  —  Bcn-j  AcOi  —  Arn-, 

f  .  -S-- — ~  tr,  +  X-4  „  .  -S  J.  „  f- 3  {  +  J_L_L1+C(n 


w 


w 


w 


w 


or,  using  (7.3.3)  and  (7.3.5) 

H  =  -  [Cn  +  hpn  -  -  n^hp  +  cw(n  -  n3)(n  -  n^)]Bp 

—  ^  [^n  +  h^n  -  0|f.  -  n-|h^  +  cw(n  -  n3)(n  - 


3M"  ~  "1 

(7.3.7) 


(7.3.8) 


.  2  . 

where  C  and  c  are  arbitrary  constants.  The  absence  of  a  term  in  f,  is 
due  to  the  orientation  of  the  parallel  sides  along  the  £  axis.  All  other 
terms  of  the  general  second  degree  polynomial  are  present  in  general.  If, 
however,  Bp  =  B^,  as  Is  true  in  the  "step  function"  option,  then  the 
quadratic  terms  vanish  and  n  is  a  linear  function  of  l,  and  n. 
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7.3.4  Variation  Between  Elements  of  a  Lifting  Strip 

The  variation  of  dipole  strength  across  the  N-Hnes,  i.e.,  the  variation 
from  one  lifting  strip  to  the  adjacent  one,  is  discussed  above.  It  remains 
to  discuss  the  variation  along  a  lifting  strip,  i.e.,  the  variation  from  one 
lifting  element  to  the  next  one  of  the  strip.  The  dipole  strength  along  the 
"top"  side  of  the  element  between  the  points  U^,  nj)  and  (s2,  n2)  (see 
figure  20)  is  obtained  by  setting  i  =  +  b^  in  (7.3.7).  The  result  is 

h(32)  =  ^(linear)  +  (Bp  -  Bg)  (cw2  +  wm32>  [J  (J-  1)]  (7.3.9) 

In  the  square  bracket  s  denotes  arc  length  along  the  side  and  L  the  total 
length  of  the  side  (L  =  d^  in  the  notation  of  seciton  7.2).  The  function 
p(linear)  is  a  linear  function  that  varies  from  the  value  of  p  at  the 
point  n-j),  which  equals  B<. .  (arc  length  of  the  point  from  the  trailing 

edge)  to  the  value  of  p  at  the  point  (e;2,  n2)  which  equals  Bp*  (arc 
length  of  the  point  from  the  trailing  edge).  On  the  adjacent  element,  the 
"bottom"  side  that  lies  between  the  points  (f;^,  n^)  and  (5.,  n-j)  is  the 
one  that  lies  along  the  side  discussed  above.  The  dipole  strength  along  this 
side  Is 

pi (41 )  =  ^(linear)  +  (Bp  -  B$)  {cw2  +  wm^ )  -  1)]  (7.3.10) 

Ignoring  the  small  gaps  between  elements  produced  by  the  projection  of  the 
input  points,  the  quantities  p(Hnear),  s,  and  L  are  Identical  In  equa¬ 
tion:  (7.3.9)  and  (7.3.10),  as  are  Bp  and  B^.  The  only  quantities  that 
ar>-  different  are  those  In  the  curly  brackets.  Here  c  and  w  correspond 
to  different  elements,  while  the  slopes  and  m^  correspond  to  dif¬ 

ferent  sides  of  different  elements.  It  is  clear  from  figure  20  that  the 
products  wm32  and  wm^  are  just  the  changes  in  the  ^-coordinate  between 
the  endpoints  of  the  side  question.  These  may  be  put  in  vector  form.  Let 
the  vector  between  the  endpoints  of  a  side  be  denoted  m.  Since  a  common 
side  of  two  adjacent  elements  Is  being  considered  (ignoring  any  higher  order 
gaps  produced  by  the  "adjustment"  process  of  section  7.2),  the  same  vector 
m  applies  to  both  elements.  Then  the  change  in  ^-coordinate  over  that  side 
is  m  '  Tp  where  as  defined  in  section  7.2,  Tp  is  the  unit  vector  along 
the  t-axis.  Now  the  change  in  dipole  strength  across  the  side  common  to 
two  elements  Is 


67 


(7.3.11) 


Ahl  =  (Bp  —  Bs){a(w2c)  +  m  •  (atE)}Cf  if-  -  1)] 

where  any  quantity  preceeded  by  A  represents  a  change  in  that  quantity.  If 
the  N-lines  are  straight  and  the  elements  are  coplanar,  aT£  =  0.  If  the 
angle  between  two  elements  is  small  Tp  Is  of  the  order  of  the  square  of  the 
angle.  Moreover,  this  angle  is  small  if  the  slope  of  the  surface  is  continu¬ 
ous  and  if  enough  elements  are  used  to  Insure  calculational  accuracy.  Thus, 
in  the  present  method  the  parameter  c  is  set  equal  to  zero  for  all  elements 
on  the  body  surface.  The  resulting  discontinuity  in  dipole  strength  between 
elements  of  an  N-line  Is  of  higher  order  than  the  other  approximations  of  the 
method  if  the  slope  of  the  body  is  continuous.  At  a  slope  discontinuity 
the  dipole  strength  can  be  made  continuous  by  having  the  N-lines  intersect 
the  line  of  discontinuity  at  right  angles,  so  that  in  •  Tp  =  0.  However, 
this  appears  to  be  generally  unnecessary  for  good  accuracy. 

One  exception  to  the  above  rule  is  the  trailing  edge.  The  local  slope 
discontinuity  is  quite  severe,  and  requiring  the  N-lines  to  be  perpendicular 
to  the  trailing  edge  Is  undesirable.  Thus,  If  only  the  on-body  dipole 
distribution  Is  considered,  there  is  a  discontinuity  of  dipole  strength  at 
the  trailing  edge,  having  the  parabolic  variation  of  the  square-bracketed 
term  In  (7.3.11)  and  thus  a  concentrated  vortex  filament  of  this  form  would 
lie  along  the  trailing  edge.  This  difficulty  is  disposed  of  by  adding  a  dipole 
term  of  the  correct  form  to  the  distribution  in  the  wake.  In  the  wake,  the 
dipole  strength  is  constant  along  N-lines  and  thus  has  the  form  of  equation 
(7.3.7)  with  Bp  =  Bs  =  0  (Ap  and  are  set  equal  to  the  actual  values 

of  Bp  and  B<.  multiplied  by  the  total  on-body  arc  lengths  of  the  respective 
N-lines).  Thus,  a  value  of  C  may  be  chosen  which  Is  proportional  to 
(Bp  —  Bj.}.  That  eliminates  the  discontinuity.  By  factoring  out  (Bp  —  B^) 

C  may  be  replaced  by  c.  In  a  manner  analogous  but  not  identical  to  the 
redefinition  involved  in  going  from  equation  (7.3.7)  to  (7.3.8).  The  resulting 
formulas  are  given  in  section  7.9,  which  deals  with  the  wake  elements. 

The  discontinuity  discussed  above,  together  with  its  remedy,  occur  only 
if  the  "piecewise  linear"  option  is  used  for  bound  vorticity.  If  the  "step 
function"  option  is  used,  then  on  an  element  Bp  —  B$  =  0,  and  the  discontinuity 
disappears.  This  is  another  simplification  connected  with  this  option. 
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7.4  Overall  Logic  of  the  Calculation  of  the  Velocity  Induced 
by  a  Lifting  Element  at  a  Point  in  Space 

The  basic  formulas  of  the  present  method  are  those  giving  the  velocity 
Induced  by  an  element  at  points  in  space.  These  are  applied  to  obtain  the 
effects  of  the  elements  at  each  other's  control  points.  For  an  element  of  a 
lifting  section  on  the  body  surface  there  are  two  kinds  of  induced  velocities, 
that  due  to  the  constant  source  density  on  the  element  and  that  due  to  the 
dipole  distribution  of  the  form  (7.3.8).  For  different  ranges  of  distance 
between  the  element  centroid  and  the  point  where  velocities  are  evaluated, 
different  sets  of  formulas  are  used.  The  three  ranges  are  denoted:  (1)  the 
far-field  or  point  singularity  regime,  (2)  the  intermediate  field  or  multipole 
regime,  and  (3)  the  near  field  or  exact  regime.  The  near-field  formulas  are 
obtained  by  an  exact  integration  over  the  elements.  Such  formulas  are  necessary 
to  obtain  the  desired  accuracy  at  points  near  the  element,  but  are  quite  time- 
consuming.  At  points  further  from  the  element  approximate  integration  formulas 
are  used  to  reduce  computing  time.  When  the  distance  between  the  element 
centroid  and  the  point  where  velocities  are  being  evaluated  exceeds  a  certain 
multiple  of  the  maximum  diagonal  of  the  element,  approximate  formulas  are  used. 
In  the  far  field,  velocities  are  calculated  directly  In  the  reference  coordinate 
system.  In  the  Intermediate  and  near  fields  the  field  point  where  velocities 
are  to  be  evaluated  is  transformed  into  the  element  coordinate  system  using 
the  transformation  matrix  (7.2.10).  Then  velocities  are  computed  in  element 
coordinates,  and  finally  the  computed  velocities  are  transformed  back  to  ref¬ 
erence  coordinates  using  the  transformation  matrix.  This  procedure  is  well 
known  and  will  not  be  discussed  further  here.  A  complete  description  is 
contained  in  reference  3. 

Now  notation  will  be  introduced  for  the  velocity  calculation.  It  is 
assumed  that  the  velocities  that  are  being  computed  are  due  to  the  j-th 
element  and  are  being  evaluated  at  the  control  point  (centroid)  of  the  i-th 
element.  Clearly,  any  point  could  be  substituted  for  the  i-th  control  point. 

The  velocity  due  to  the  constant  source  density  is  denoted  ?. ..  In  element 

*  J 

coordinates  It  has  components  V  (source),  V  (source),  and  V  (source), 

X  jr  c. 

1  .e., 

=  Vx(source)  +  V^(source)  +  Vz(source)  (7.4.1) 
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For  a  general  quadrilateral  element  of  a  nonlifting  section,  this  source 
velocity  is  the  only  induced  velocity,  and  it  is  computed  by  the  formulas  of 
reference  3  in  all  three  ranges.  For  a  trapezoidal  element  of  a  lifting 
section  the  calculation  of  source  velocities  is  the  same  as  for  a  nonlifting 
element  in  the  far  field  and  the  Intermediate  field  (a  trivial  difference  is 
the  use  of  normalized  area  moments).  In  the  near  field  advantage  can  be 
taken  of  the  fact  that  the  element  is  a  trapezoid  to  shorten  the  formulas 
and  conserve  computing  time. 

To  develop  formulas  for  the  velocity  induced  by  the  dipole  distribution 
on  an  element,  some  additional  notation  is  required.  Furthermore,  it  simpli¬ 
fies  the  development  to  consider  the  velocity  potential  initially  rather 
than  the  velocity  components.  The  potential  due  to  the  dipole  distribution 
on  the  element  at  points  of  space  is  obtained  by  integrating  over  the  element. 
Namely, 

t  =  JT ({’(dipole)  nU,n)  d£dn  (7.4.2) 

E 

where  tiU.n)  Is  given  by  (7.3.8), where  the  integration  is  over  the  area  of 
the  element,  and  where  <j>(di pole)  is  the  potential  of  a  unit  point  dipole 
with  axis  normal  to  the  element,  i.e., 

<{.(dipole)  =  —  ■  - gyTj-  (7.4.3) 

C(z  -c)2  +  (y  -n)2  +  Z2] 

Here  (x,  y,  z)  Is  the  point  where  the  potential  and  velocity  afe  being 
evaluated  expressed  In  the  coordinate  system  of  the  element.  Now  define 
the  auxiliary  potentials 

*pq  "If  ^^diP°1e^Pr'q  d^dn  (7.4.4) 

E 

where  p  and  q  are  integers.  Now  from  (7.3.8),  (7.4.2),  and  (7.4.4),  the 
potential  of  the  element  is 

$  *  m pBp  —  (7.4.5) 
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where 


=  w  ^11  +  Voi  “  n3^10  “  n3hF^OO  +  CW*c] 

=  w  ^11  +  Voi  “  n1  ^10  ~  r|lhS<J>00  +  cw»c]  (7.4.6) 

$c  =  $02  ~  ^nl  +  n3^01  +  nln3l),00 

As  stated  in  the  previous  section,  the  term  $  is  not  currently  used  for 

lifting  elements  on  the  body.  For  completeness,  it  is  included  in  the  formu¬ 

lation,  and  equations  are  qiven  in  the  subsequent  sections.  These  last  are 
needed  for  wake  elements  in  any  event.  The  velocity  due  to  the  dipole  distri¬ 
bution  is 


V^j(dipole)  =  -  Vtf,  (7.4,7) 

where  v  denotes  the  gradient  operator.  In  element  coordinates  this  is 


From  (7.4.5)  and  (7.4.7)  it  is  clear  that 

V, ^dipole)  -  Bp  V<$>  *  Bs 


where 


W  ‘  "%* 


-  »7*s 


(7.4.8) 


(7.4.9) 


(7.4.10) 


The  desired  velocities  are  these  .,  \[V ,  and  In  the  far  field 

these  are  calculated  directly.  In  the  near  and  Intermediate  field  the  source 
velocity  is  evaluated  directly,  but  the  dipole  velocities  are  broken  up  into 
separate  terms  in  the  manner  of  (7.4.6).  Thus  to  evaluate  the  dipole  velocities, 
formulas  are  needed  for  the  derivatives  of  $qq,  4>-jq,  $qj  ,  and  (j^. 

These  formulas  are  presented  in  the  following  sections. 


As  mentioned  above,  the  Integrals  (7.4.4)  can  be  evaluated  analytically 
and  the  resulting  expressions  differentiated.  This  is  what  is  done  for  the 
near  field  (section  7.7).  The  resulting  expressions  are  quite  Involved  and 
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A 


time-consuming  to  evaluate.  To  save  computing  time,  approximate  formulas  are 
used  when  the  field  point  is  some  distance  from  the  element.  This  is  accom¬ 
plished  by  means  of  a  multipole  expansion.  The  basic  idea  is  to  expand 
<f> (di pole)  from  (7.4.3)  in  a  two-dimensional  Taylor  Series  about  £  =  n  =  0. 

This  process  is  a  standard  development  in  the  textbooks.  The  result  is 

2 

<+. ( di pol e)  =  FQg(x,y,z)  F-jQ(x,y,z)£  +  FQ^(x,y,z)n  +  Fgg(x,y,z)£ 

t  Fn'x,y,2)tn  ♦  F02(x,y,2)n  *  ^  ^  FJx'y-zk  "  *  — 

n  m 

(7.4.11) 

where  the  F  are  the  derivatives  of  <*>(dipole)  at  the  origin  of  element 
nm 

coordinates  and  are  independent  of  £  and  n.  When  (7.4,11)  is  inserted  into 
(7.4.4),  the  Fnm(x,y,z)  are  taken  out  of  the  integral,  the  remaining 
integrals  are  of  the  form  (7.2.24)  and  are  thus  moments  of  the  area  of  the 
element,  which  can  be  normalized  by  division  by  the  appropriate  powers  of  t. 

In  the  intermediate  field  the  expansion  (7.4.11)  is  used  through  the 
second-order  terms,  F2Q,  F^,  Fq£.  In  the  far  field  only  the  initial,  zero 
order,  term  is  used.  It  is  clear  from  the  form  of  (7.4.11)  that  Fqq  is  the 
potential  of  a  unit  point  dipole  at  the  origin  of  element  coordinates  (centroid). 
In  the  far  field  every  auxiliary  potential  (7.4.4)  is  a  multiple  of  the  point 
dipole  potential  and  thus  so  are  the  combined  potentials  (7,4.6).  Thus, 

Induced  velocities  In  the  far  field  may  be  expressed  directly  in  reference 
coordinates  using  the  well-known  formulas  for  a  point  dipole. 

The  above  discussed  only  the  dipole  velocity,  but  the  same  procedure  is 
followed  for  the  source  velocity.  In  fact,  the  development  for  this  case  is 
given  in  detail  in  reference  3. 

7.5  Par-Field  Formulas  for  the  Velocity 
Induced  by  a  Lifting  Element 

First  calculate  the  distance  rQ  between  the  centroid  of  the  element 
and  the  field  point  where  velocities  are  to  be  calculated.  If  the  reference 
coordinates  of  the  centroid  are  xQ,  yQ,  zQ  and  the  reference  coordinates  of 
the  field  point  are  x',  y',  z',  then 
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r0  =  ^  (*'  “  V  +  *y'  "  yo'  +  (z‘  “  V  ( 

Now  test  r0/t,  where  t  Is  the  maximum  diagonal  of  the  element.  If 

rQ/t  >  p,  ( 


(7.5.1) 


(7.5.2) 


where  p-|  is  a  certain  criterion,  then  the  far-field  formulas  are  used. 
Currently  is  set  equal  to  4.0.  The  far-field  formulas  calculate 
velocities  directly  in  reference  coordinates.  First  define  the  vector 

r0  =  (x*  -  xQ)t  +  (y*  -  yQ)T  +  (z‘  -  z0)TT  (7.5.3) 


The  source  velocity  is 


The  dipole  velocities  are 


where 


rr  .  t  *00  ro 

1j '  “2T  r° 

o 


T<0  -  -QFt> 


V'.V  =  Q«-ff 

ij  S 


~  w  ”T  *11  n3hLI00  +  CW^  *02  +  ’’’lVoO^ 
o 


QS  =  w  *11  nlhRI00  +  cw^t2l02  +  nln3I00^ 

'  o 


(7.5.4) 


(7.5.5) 


(7.5.6) 


and  where 


/n  •  r  \  r 

-  3M 


(7.5.7) 


It  will  be  recalled  that  rT  is  the  unit  normal  vector  of  the  element  (n  is 

not  connected  w<th  the  field  point)  and  that  I  denotes  normalized  moments 

nm 

as  given  by  (7.2.24).  A  comparison  of  (7.5.6)  with  (7.4.6)  shows  that  the 
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4qi  and  4> ^ Q  terms  have  dropped  out  because  they  are  multiplied  by  the  zero 
value  moments  I q -j  and  I^. 

7.6  Intermediate-Field  or  Multipole  Formulas  for  the 
Velocity  Induced  by  a  Lifting  Element 

If  ro/t  <  transform  the  reference  coordinates  x ' ,  y',  z'  of  the 
field  point  by  the  transformation  matrix  to  obtain  element  coordinates  x,  y,  z 
of  the  field  point.  Now  perform  another  test.  If 

rQ/t  >  p2  (7.6.1) 

where  p2  is  another  input  criterion,  which  is  currently  taken  as  p2  =  2.5, 
then  the  multipole  formulas  are  used.  The  dipole  velocities  are  taken  in  the 
form  (7.4.10),  which  means  that  derivatives  of  all  quantities  in  (7.4.6)  must 
be  calculated. 


First  define  direction  cosines 


a  =  — , 
r 

0 

e 

=  iL 

r  ’ 

0 

z 

7  =  _. 

o 

(7.6.2) 

Next  define  certain 

"derivative 

functions"  as 

follows: 

First  Order 

u  =  —  a, 

X 

U  “ 

y 

-  3, 

uz  =  -  r 

(7.6.3) 

Second  Order 

u  =  3a^  -  l, 

XX 

u  = 
xy 

3aP, 

u  =  3P2  ~  1 

2 

(7.6.4) 

u  =  3a  y 
xz 

u  = 

yz 

337, 

uzz  '57-1 

Third  Order 

<W  '  *'5  - 

u  - 

rxy 

330  -  ^), 

U  -  37(1  -  3a2', 

xxz 

">U  '  *<l  -  5f2) 

u 

xyz 

-  15apy 

u  =  3a(l  -  372) 

(7.6.5) 

uyyv  ’  3MJ  - 

u  = 

yyz 

3rd  -  532) 

v  =  }e<l  -  ^ 
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Fourth  Order 


XXXX 


=  9  -  90a2  +  105a4 


\xxy  *  “  5) 


uxxxz  =  -  3) 

Uxxyy  =  3  ~  ^(o2  +  P2)  +  105C^&2 


Uxxyz  =  W7(T0?  -  1) 


=  3  -  15(0^  +  72)  +  105a2 y2 


XXZZ 


>W  *  -  3) 

>W  ’  w(7fs2  ~  « 


\yzz  '  *  1) 

2 

u  =  9  -  900  +  1050 

yyyy  ^  ^ 


W  -  1^<TT  -  3) 

“jry^Z  "  3  -  15(P2  ♦  32)  +  lOSe2?2 


Then  the  source  velocity  components  are 


V,  (source)  = 

r 

o 


V  (source)  =  -77 
y  2 

J  r 

o 


Vz( source)  - 

r 

o 


(7,6*6) 


I00Ux 


_u  +  21.  ,u  +  I  u  ] 
0  xxx  11  xxy  02  xyy 


KO  113 

“  I00Uy  ”  2  (r^)  [l20Uxxy  +  2IllUxyy  H'  I02Uyyy‘1 


(7.6.7) 


I00Uz  2  (r0)  ^*2 


artiVi*- 

+  21,  .u  "  ■  +  I  'a  J 
0  xxz  11  xyz  02  yyz 


These  are  identical  to  the  miltipole  formulas  of  reference  3  with  a  slightly 
different  notation. 


Specific  formulas  for  the  derivatives  cf  the  various  dipole  potentials 
<t>pq  appearing  in  (7.4.6)  are  given  on  the  following  page.  To  illustrate  the 
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Dipole  Derivatives  for  Multipole  Expansion 


Far -Fie  id  1st  Order 


2nd  Order 


■  -  4  V«  (¥){KHA  *  W]  *<kf[  h?: 


+  21,,  u  +  I_,j 
xxxz  11  xacy i  02 


»  l! 

:  xyyzj  j 


-  $  j:ooV  -  (tyiKKL  t  *  (y  [  WW  *  *  Vjm]  !  <7  -6-8> 

o 

■  ^3  |X00UZZ  ”  t  ^0 (u^]  +  (“)  [  I20Uxxzz  +  2IllUxyzz  +  ^“yyzz]  | 

o 

-  ^3  ~  (  r^)  [  Wxxz  +  ^^xyz]  +  (r^)  [  I30Uxxxz  +  2I21uxxyz  +  I12uxyyz]  ) 

o 

~  ^3  "(r^)  [X20Uxy7.  +  IllUyyz]  +  (r^)  [  ^“xxyz  +  2I21uxyyz  +  I12uyyyz]  ]  ^ 

O 

“  ^3  jn.^2  “  (r~)  [*20uxzz  +  ^lluyzz]  ’  (r^)  [X3Cuxxzz  +  2I21uxyzz  +  ^"yyzz]  ] 

o 

"  ^3  {*oKz  _  (r^)  [  Ill“xxz  +  I02Uxyz]  +  (l~)  [^^xxxz  +  2I12uxxyz  +  I03Uxyyz]  | 

:  “  3  ~  (^)  [^xyz  +  I02Uyyz]  +  (^)  [-^Axyz  +  2I12Uxyyz  +  ^yyyJ  j  t7*6-10^ 

1  "  3  | ^oinL  ~  (r^)  [IUuxzz  +  I02Uyzz]  +  (r^)  [^l^cxzz  +  2I12uxyzz  +  X03uyyzz]  | 

:  "  ^3  |IllUxz  “  (?^)  [I21Uxxz  +  I12Uxyz]  +  (i")  [SlUxxxz  +  2I22uxxyz  +  I13Uxyyz]  | 

;  "  ^3  KlV  “  (^)  [^“xyz  +  *12uyyz]  +  (?^)  [X31uxxyz  +  2I22uxyyz  +  X13uyyy*]  |  t7*6*11) 

ro 

ll  |  2  *1  ) 

-  ^3  |Xlluzz  (rQ)  [^^xzz  +  XL2uyzz]  +  (rQ)  [X31Uxxzz  +  2X22uxyzz  +  I13UyyzzJ  j 
ro 

=  3  X02Uxz  ~  (~)  [ri2Uxxz  +  X03Uxyz]  (  r0 )  [I22Ux<xz  +  2X13Uxxyz  +  X04uxyyz]  | 

ro 

=  -  j  I  U  ~(— )  T  I,0U  +  I„.U  ]  +  (~)  f  I„0U  +  2I..U  +  I  ,  U  l!  (7.6.12) 

3  (  02  yz  \ rQ  /  [  12  xyz  03  yyzj  Vr0/  I  22  wz  13  xyyz  04  yyyzj  j 
ro 

“  “  3  |102uzz  “  (r^)  [I12uxzz  +  X03Uyzz]  +  (r^)  [*22Uxxzz  +  2I13Uxyzz  +  X04uyyzz]  j 
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development, the  terms  containing  the  moments  I^g  and  are  written  down. 

They  are  then  crossed  out  to  show  that  such  terms  need  not  be  calculated 
because  =  Ig-j  =  0.  (Inclusion  of  these  terms  generalizes  the  formulas 
to  the  case  where  the  centroid  is  not  used  as  origin.) 

7.7  Near-Field  Formulas  for  the  Velocity  Induced 
by  a  Lifting  Clement 

If  r  /t  <  p£»  the  near-field  or  exact  formulas  are  used  to  compute 
Induced  velocities.  The  calculation  starts  with  the  element  coordinates 
x,  y,  z  of  the  field  point  and  the  geometric  quantities  associated  with  the 
element  that  are  discussed  in  section  7.2.  In  particular,  the  corner  point 
coordinates  nk»  k  =  1,  2,  3,  4  are  needed,  together  with  the  width  w 
from  (7.2,14),  the  slopes  rn^  and  m^  from  (7.2.15),  the  intercepts  b32 
and  b^-j  from  (7,2.20),  the  maximum  diagonal  t  from  (7.2.21),  and  the  side 
lengths  d^i  d^t  d34,  d^, ,  from  (7.2.22).  These  quantities  are  illustrated 
in  figure  20.  In  addition  to  the  basic  near-field  equation,  there  are  special 
limiting  formulas  for  small  values  of  rQ  and  z.  However,  the  basic  near- 
field  formulas  are  used  in  the  large  majority  of  cases. 
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and 


I) 

run 


(r  +  r  Y~  -  (f 

ill  ll  1 


(7.7.6) 


ran  -  12,  2'j ,  3**,  4l 

Now  In  terms  of  the  above  functions  the  source  velocities  and  dipole  potential 
derivatives  needed  1 n  (7.4.6)  can  be  written. 


78 


The  source  velocities  are 


V  (source)  -  — 
x 


L<32)  + 


1  +  ra 


32 


V1  +  m! 


■I - L(^l) 


2 

*1 


V  (source) 

y  ' 


L(“)  +  ipM  +  mX  l(52)  _  V  L(U) 


V 


1  +  m 


32 


y 


1  +  m. 


(7.7.7) 


‘4l 


V (source) 


(32)  (32)  (41)  (41) 

2  *3  1  “  T4 


To  evaluate  the  dipole  potential  derivatives,  the  derivatives  of  Vx  and 

Vy  are  needed  (since  Vz  *  Its  derivatives  are  exactly  a  potential 

derivative).  The  derivatives  of  V  and  V  are 

x  y 


dV  (source) 

1 

+ . 

1 

&x 

V1  +  m 32 

3x 

\ 

2 

m4l 

c)x 

SVx( source) 

1 

aiP2>  t 

_ 1_ 

— - ^ - 

3y 

oy 

V1  +  m32 

\l~ 

"  2  ' 

oy 

dV  (source) 

X 

1 

ai,<32> , 

1 

a/>l> 

OZ 

/ - 2~ 

V1  +  m32 

V1  +  m4] 

02 

(7.7.8) 

<)v  (source) 

-Xy - -  = 

Sl(12)  ,  *P4> 

+ 

Sl,<32> 

r,14i 

*,<4l> 

dx 

c: .  3x 

\A*4 

3x 

V1  +  4 

dx 

SV  (source) 

-JL _ -  = 

aLl12'  .  5L<31l> 

,  "52 

9L<32> 

“41 

»L<4l> 

<3y  3y 

dy 

£y 

V1  +  ra32 

\A  *  4 

(source) 

y 

—  ......  s. 

.  SiP-l 

,  rc32 

ai(32> 

m4l 

c*L<Ul> 

<3z 

02  02 

VL  +  m32 

32 

V1 + 4 

dz 
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Now  the  potential  derivatives  are  as  follows. 


?*00  . 
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OX 

3x 
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dy 

Sy 
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y 
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S*01 
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- 
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5*01 
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^10 
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oy 
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(7.7.9) 


(7.7.10) 


(7.7.11) 


Evaluation  of  the  derivatives  of  and  requires  certain  auxiliary 

quantities  and  and  their  derivatives.  Thus  define 
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7.8  Some  Alternate  Near-Field  Formulas  for  Use  in  the 

Plane  of  the  Element 

If  the  point  where  the  velocity  induced  by  an  element  is  being  calculated 
lies  In  the  plane  of  the  element,  i.e.,  If  2=0,  there  may  be  numerical 
difficulties  In  the  evaluation  of  the  formulas  of  section  7.7  for 
^OO  =  ^z ( source)  and  Its  2-derivatlve.  To  avoid  possible  difficulty  special 
formulas  have  been  derived  for  this  case. 


If 


| z/t 1  <  0.001  (7.8.1) 

the  point  (x,y,z)  is  considered  to  be  in  the  plane  of  the  element  and  z  is 
set  equal  to  zero-  Now  Vz(source)  is  2n  for  points  Inside  the  element  and 
zero  for  points  outside.  Some  tests  for  this  condition  have  encountered 
problems  of  numerical  significance.  The  following  tests  are  currently  used. 
First  define 


h32  =  m32  n3^  C3^ 

^4i  =  m4i  (y  — ■  x  ^ i ) 


(7.8.2) 


Then  a  point  is  Inside  the  element  if  all  three  of  the  following  tests  are 
satisfied  and  outside  if  any  one  is  not  satisfied. 


rQ/t  <  1/2 

h32h41  <  0  (7.8.3) 

(y  —  n-j )  (y  —  03)  <  0 

The  velocity  V2(source)  is  simply  set  equal  to  2v  or  to  zero  as  appropriate. 


Numerical  difficulty  can  be  encountered  in  the  evaluation  of  the 
z-derlvative  of  4>q0  when  the  point  (x,y,0)  is  on  an  extension  of  a  side  of 
an  element.  This  condition  can  be  determined  by  testing  the  above-defined  h's 
and  the  y  —  n.  Specifically,  the  point  is  considered  to  be  on  a  side  if  any 
of  the  following  tests  are  satisfied  (refer  to  figure  20  for  element  geometry): 
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(7.8.4) 


Point  on  Side  12  if  |y  —  |/t  <  0.001 
Point  on  Side  23  if  ! I /t  <  0.001 
Point  on  Side  34  if  |y  -  n^l/t  <  0.001 
Point  on  Side  41  if  |h^^|/t  <  0.001 


If  none  of  the  above  tests  are  satisfied,  the  formulas  of  section  7.7  are 
used  for  the  z-derivative  of  4>qq.  Only  one  of  the  conditions  (7.8.4)  can  be 
true.  If  this  occurs,  then  the  following  substitution  is  made  in  equation 
(7.7.9). 


Side  12: 
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Side  23: 


Side  34: 
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Side  41 : 
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The  remaining  two  T  derivatives  of  equation  (7.7.9)  are  evaluated  by  the 
formulas  of  section  7.7. 


7.9  The  Velocity  Induced  by  a  V/ake  Element 


In  the  wake  the  dipole  strength  Is  constant  along  N-lines,  as  illustrated 
In  figure  21.  The  form  of  the  dipole  distribution  on  a  wake  element  is 
obtained  by  setting  Bp  =  =  G  in  equation  (7.3.7).  Specifically, 

p.  =  ^[(Ap  —  A<.)n  +  A^n-j  —  Apn-j]  +  C(n  —  n3)(n  —  ii )  (7.9.1) 

Denote  by  L  (total)  the  total  arc  length  along  an  N-line  from  the  trailing 
edge  around  the  section  curve  of  the  body  and  back  to  the  trailing  edge.  This 
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arc  length  Is  computed  in  a  manner  similar  to  (7.2.23),  namely 


Lp  (total)  =  ^dp 
Ls  (total)  = 


(7.9.2; 


where  the  suns  are  over  all  on-body  lifting  elements  of  the  strip  lying 
between  the  two  N-lines.  Now  from  the  form  of  the  dipole  distribution  shown 
in  figure  21  it  is  clear  that  the  constant  values  Ap  and  A<.  assumed  along 
the  N-llnes  in  the  wake  are  equal  to 


Ap  =  BpLp  (total ) 

=  BsLs  (total) 

Thus,  velocity  potential  due  to  a  wake  element  has  the  form 


(7.9.3) 


$  =  4>pBp  (7.9.4) 

just  as  in  equation  (7.4.5).  Here,  however. 


*F  =  5  [*0!  -  Voo’jLF  (total>  +  «c 

't’S  _  ^  *^1  ;OO^S  ^ total)  +  (7,9.5) 


^c  =  ^02  ^nl+  n3^01  +  nln34>00 

These  replace  (7.4.6).  The  dipole  velocity  is  given  as  before  (see  (7.4.7), 
(7.4.9),  and  (7.4.10)  by 


(7-9-6) 

where 

"V  ’ll’ s  (7-9-7> 

To  evaluate  (7.9.6)  in  the  near  and  intermediate  field,  the  derivatives  of 
‘*’02’  ^OT  and  ^00  are  evaluated  by  the  formulas  of  sections  7.7  and  7.8. 
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In  the  far  field  the  formulas  for  the  dipole  velocities  due  to  a  wake 
element  are 


where 


QF  =  w 


t2  Too 

7~ 

0 


n^Lp(tot) 


and  where  as  before  (see  (7.5.7)) 


(7.9.8) 
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(7.9.9) 


(7.9.10) 


There  is  no  source  density  on  wake  elements  and  no  source  velocities  are 
computed . 


As  discussed  in  section  7.3.4,  the  values  of  c  on  wake  elements  are 

not  zero  if  the  "piecewise  linear"  option  for  bound  vorticity  is  usee. 

Instead,  the  value  of  c  on  the  first  wake  element  is  determined  to  avoid 

a  discontinuity  in  dipole  strength  at  the  trailing  edge.  Values  of  c  on 

the  remaining  wake  elements  are  chosen  to  eliminate  discontinuities  between 

adjacent  wake  elements  along  the  lifting  strip.  Let  superscript  (1)  denote 

quantities  associated  with  the  first  on-body  element  of  a  lifting  strip  and 

superscript  u  denote  quantities  associated  with  the  last  on-body  element  of 

the  strip.  Similarly,  the  superscripts  wl,  w2,  etc.  denote  the  first  wake 

element,  second  wake  element,  etc.  of  the  same  lifting  strip.  The  important 
{ w  1  ) 

value  of  c  is  c  ,  i.e.,  the  one  for  the  first  wake  element.  It  is 

computed  from 


.(wl)  _ 


„(u)[w(u)c(u) 


m3?^ 


„<V»c<n 


(7.9.11) 


where  the  quantities  w,  m^.  have  their  usual  meaning.  Values  of  c 

fer  the  remaining  wake  elements  are  obtained  from 

c<“' >[„(»' )]2  .  c(w2)rw(w2)12  ,  c(«3)[„(w3)]2  „  .  .  . 
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7.10  Option  for  a  Semi-Infinite  Last  Wake  Element 


In  most  cases  of  Interest  the  trailing  vortex  wa  *  extends  to  infinity. 

To  facilitate  accounting  for  this  condition,  provision  has  been  made  for  con¬ 
sidering  the  last  element  of  the  wake  to  be  semi-infini fe.  A  finite  element  of 
the  sort  shown  in  figure  20  Is  formed  at  the  end  of  the  wake,  Including  all 
the  geometric  quantities  of  section  7.2.  The  induced  velocity  calculation  for 
this  element  Is  performed  using  the  origin  of  coordinates  appropriate  to  the 
finite  element,  but  the  formulas  used  to  calculate  induced  velocities  are 
appropriate  for  the  semi-infinite  element.  Naturally,  all  points  in  space  are 
in  the  "near  field"  with  respect  to  a  semi-infinite  element,  so  it  is  the 
formulas  of  section  7.7  that  apply.  These  'ormulas  are  modified  by  setting 


m3?  =  0 


{7.10.1} 


2  3 

This  yields  immediately 

u1»  y  i  •  -l4*  ^4*  uncha nged  (7.10.2) 

S3  ~f'2  =  y3  =  >2  =  0 
The  log  functions  (7.7.3)  and  their  derivatives  (7.7.6)  are  replaced  by 


03  —  ~  *  1 1 


(41) 


unchanged,  all  derivatives  unchanged  (7.10.3) 


=  0,  all  derivatives  equal 


zero 


(12)  (34)  r4  ?’4^ 

-L(  )  +  l(  *  =  l0s  rp-or-xpr 


(7.10.4) 
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small  compared  to  the  source-velocity  matrices.  Each  of  the  velocities 
(7.11.1)  represents  the  velocity  due  to  a  dipole  distribution  on  the  strip 
that  Is  unity  on  one  N-llne  and  zero  on  the  other  with  a  linear  "spanwise" 
variation  in  between. 


The  characteristic  onset  flow  velocities  due  to  a  strip  are 


=7<s)  +  v<F) 

vik  vik  vik 
vik  2  LVik  vik  J 


(7.11.2) 


The  first  velocity  of  (7.11.2)  is  that  due  a  dipole  distribution  on  the  strip 

that  is  constant  in  the  "spanwise"  direction.  The  second  velocity  is  that  due 

to  a  dipole  distribution  that  varies  llnearlly  in  the  "spanwise"  direction  and 

has  zero  value  at  "midspan  ",  These  velocities  are  used  to  form  the  basic 

kl 

circulatory  onset  flows  V:  . 


If  the  "step  function"  option  for  bound  vorticity  is  used  (section  6.3), 
the  proper  form  of  the  dipole  distribution  is  simply  constant  in  the  "spanwise" 

■Ao) 

direction  over  a  lifting  strip,  and  the  velocity  V-k'  is  precisely  the  onset 
flow.  Thus,  for  this  option,  the  vorticity  onset  flows  are 

vjk)  =  ,  k  =  1,  2 . L  (7.11.3) 

The  above  yields  L  onset  flows,  each  of  which  corresponds  to  a  unit  value 
of  the  "streamwise"  dipole  derivative  B  on  one  lifting  strip  and  zero  values 
of  B  on  all  other  lifting  strips.  (Recall  that  the  "streamwise"  derivative 
of  dipole  strength  is  essentially  the  value  of  the  bound  vorticity.)  No 
special  handling  is  required  at  the  ends  of  the  lifting  section. 


The  machinery  for  the  "piecewise  linear"  option  for  bound  vorticity  is 
somewhat  more  complicated.  The  "spanwise"  variation  of  the  "streamwise" 
dipole  derivative  B  (bound  vorticity)  is  linear  in  the  "spanwise"  direction 
across  a  lifting  strip.  Thus,  t' e  velocity  at  the  i-th  point  (control  point 
or  off-body  point)  due  to  the  k-th  strip  Is 

*, (swp  k)  .?<;>  \  (7.H.4) 
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where  wk  is  the  " 
derivative  of  B, 
k-th  lifting  strip, 
through  Bk_ -j ,  Bk, 


spanwise"  width  of  the  strip,  B'  is  the  "spanwise" 
and  subscripts  k  denote  quantities  associated  with  the 


The  derivative 


Bk  is  evaluated  by  a  parabolic  fit 
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(7  11.5) 


'k "  wk  +  +  wk+]) 


wk  4  wk-l 
lwk  4  wk+l 


Then  (7.11,4)  is  approximated  numerically  by 

V-  (strip  k)  -  V<»)8k  ♦  V<J>  [0kBk.,  ♦  EkBk  t  FkBw]  (7.11.6) 


The  velocity  (7.11.6)  contains  values  of  the  "streamwise"  dipole  derivative  B 
for  three  consecutive  strips.  However,  a  proper  circulatory  onset  flow  is 
proportional  to  the  value  of  B  on  a  single  strip.  Since  each  Bk  enters 
(strip  k)  for  three  consecutive  strips,  its  three  contributions  may  be 
summed  to  give  the  basic  vorticity  onfet  flow. 
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(7.11.7) 


In  performing  the  above  parabolic  fit  (7.1.. 6),  the  values  of  the  function 
B  to  be  fit  are  of  course  the  values  of  bound  vorticity  on  the  strips.  Each 
of  these  has  been  associated  with  an  abscissa  or  "Independent  variable"  that 
expresses  the  spanwise  position  of  each  strip.  Differences  of  these  abscissas 
appear  as  combinations  of  the  widths  wk>  Calculation  of  the  wk  is  not 
obvious,  because  in  general  the  "span"  or  width  of  each  strip  is  not  constant 
but  varies  in  the  "chordwise"  direction.  Accordingly,  it  was  decided  to  input 
the  quantities  necessary  to  deduce  the  spanwise  positions  of  the  lifting 
strips.  The  input  quantities  consist  of  a  set  of  widths  wk  for  all  lifting 
strips.  If  a  strip  is  truly  of  constant  width,  it  is  natural  tc  input  that 
width.  If  the  strip  varies  In  width,  some  average  value  must  be  input  as 
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where  is  the  "spanwise"  width  of  tl> .  strip,  R‘  is  the  "spanwise" 
derivative  of  B,  and  subscripts  k  denote  quantities  associated  with  the 
k-th  lifting  strip.  The  derivative  is  evaluated  by  a  parabolic  fit 
through  , ,  B^,  and  Vr  Specifically,  define 
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(7  11.5) 


Then  (7.11.4)  is  approximated  numerically  by 

h  k>  *  ^k)Bk  +  tVk-1  +  Ekek +  FkVd  (?-n.6) 

The  velocity  (7.11.6)  contains  values  of  the  "streamwise"  dipole  derivative  B 
for  three  consecutive  strips.  However,  a  proper  circulatory  onset  flow  is 
proportional  to  the  value  of  B  cn  a  pi, ogle  strip.  Since  each  enters 
IT-  (strip  k)  for  three  consecutive  strips,  its  three  contributions  may  be 
summed  to  give  the  basic  vorticity  onset  flow. 
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In  performing  the  above  parabolic  fit  (7.11.6),  the  values  of  the  function 
B  to  be  fit  are  of  course  the  values  of  bound  vorticity  on  the  strips.  Each 
of  these  has  been  associated  with  an  abscissa  or  "independent  variable"  that 
expresses  the  spanwise  position  of  each  strip.  Differences  of  these  abscissas 
appear  as  combinations  of  the  widths  w^.  Calculation  of  the  is  not 
obvious,  because  in  general  the  "span"  or  width  of  each  strip  is  not  constant 
but  varies  in  the  "chordwise"  direction.  Accordingly,  It  was  decided  to  input 
the  quantities  necessary  to  deduce  the  spanwise  positions  of  the  lifting 
strips.  The  input  quantities  consist  of  a  set  of  widths  w^  for  all  lifting 
strips.  If  a  strip  is  truly  of  constant  width,  it  is  natural  to  input  that 
width.  If  the  strip  varies  in  width,  some  average  value  must  be  inout  as 
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the  for  that  strip,  and  this  average  is  decided  upon  by  the  user. 

These  are  used  only  in  performing  the  parabolic  fit.  To  facilitate 
fitting  at  the  first  and  last  strips  of  a  lifting  section,  it  was  decided 
originally  to  input  widths  for  ficticious  strips  adjacent  to  the  first  and 
last  strips  of  the  section.  Thus,  if  the  strips  of  a  lifting  section  were 
input  from  left  to  right,  the  table  of  would  consist  of  the  following 
sequence:  a  value  of  for  a  ficticious  strip  to  the  left  of  the  first 
lifting  strip,  the  values  of  for  the  lifting  strips  of  the  section  in 
order  from  left  to  right,  and  finally  a  value  of  w^  for  a  ficticious  strip 
to  the  right  of  the  last  strip  of  the  section.  Thus,  if  the  section  has  L 
lifting  strips,  L  +  2  values  of  w^  are  input.  This  is  still  the  format 
of  the  input.  However,  for  certain  frequently-occurri nq  situations,  the 
program  overrides  the  input  and  puts  in  a  predetermi ned  value  of  w^.  In 
fact,  it  is  only  for  the  "extra  strip"  condition  described  below  that  input 
values  of  w^  corresponding  to  ficticious  strips  are  actually  used  in  th^ 
calculations . 

Physically,  a  lifting  section  may  end  in  various  ways,  some  of  which 
involve  logical  difficulties  in  the  basic  potential-flow  model  (section  5.3). 
The  various  ways  a  lifting  section  may  end  require  various  procedures  for 
performing  the  parabolic  fit  of  the  piecewise-1 i near  vorticity  option.  These 
procedures  are  outlined  below.  In  future  work  perhaps  still  other  procedures 
will  be  required  for  situations  that  are  unanticipated  at  present. 

Sometimes  a  single  lifting  portion  of  a  three-dimensional  configuration 
is  divided  into  two  or  more  lifting  sections.  This  may  be  done  to  concentrate 
elements  in  a  certain  region,  as  shown  in  figure  23,  or  it  may  be  done  simply 
for  convenience.  In  this  case  the  division  into  sections  is  purely  logical 
rather  than  physical,  and  the  bound  vorticity  distribution  should  vary  smoothly 
from  one  section  to  another.  As  regards  the  parabolic  fit,  the  last  lifting 
strip  of  the  first  section  and  first  lifting  strip  of  the  second  should  be 
regarded  as  adjacent  strips  of  a  single  lifting  section  and  the  fit  performed 
accordingly.  This  situation  has  been  designated  "continue"  in  the  method. 

!f  the  lifting  section  has  a  physical  ending  in  the  fluid,  such  as  a 
wing  tip,  the  bound  vorticity  strength  must  fall  to  zero.  A  ficticious  logical 
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SECTION  2 


COMMON  N  -  LINE 


SECTION  I 


Figure  23.  An  example  of  division  of  a  single  physical  lifting  portion 
of  a  body  into  two  lifting  sections. 

strip  is  imagined  adjacent  to  the  tip  In  the  fluid  (figure  24a).  The  bound 
vorticlty  slope  at  the  midspan  of  the  strip  of  elements  adjacent  to  the  tip 
(figure  24)  Is  obtained  by  fitting  a  parabola  through  the  value  on  the  strip 
itself,  the  value  on  the  next  strip  inboard,  and  a  zero  value  at  the  midspan 
of  the  ficticious  logical  strip.  Various  assumptions  about  t.he  width  of  the 
ficticious  logical  strip  were  tried,  and  it  was  concluded  that  taking  its 
width  equal  to  that  of  the  lifting  strip  adjacent  to  the  tip  is  about  as  good 
a  choice  as  any,  and  this  has  been  built  into  the  program  as  an  override  to  any 
Input  value.  A  zero  width  of  the  ficticious  strip  has  a  certain  appeal,  because 
in  the  limit  of  infinite  element  number  the  bound  vorticity  must  be  zero  right 
at  the  tip.  However,  this  choice  leads  to  poor  results.  This  type  of  end  to 
a  liftlnq  section  is  denoted  "normal." 

If  a  lifting  section  ends  di  u  positive  symmetry  plane  of  the  flow 
(figure  24b),  the  proper  procedure  is  obvious.  Physically,  there  is  a  strip 
adjacent  to  the  last  strip  of  the  section  on  the  other  side  of  the  symmetry 
plane,  and  these  two  strips  have  equal  widths  and  equal  values  of  bound 
vorticity.  The  parabolic  fit  is  performed  accordingly. 
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Figure  24.  Special  procedures  at  the  ends  of  a  lifting  section  for  the 

parabolic  fit  used  with  the  piecewise  linear  vorticity  option, 
(a)  Wing  tip.  (b)  Positive  symnetry  plane. 

If  there  is  an  extra  strip  of  elements  adjacent  to  the  end  of  a  lifting 
section,  as  described  in  section  6.4,  the  width  of  the  extra  strip  is  input 
as  the  last  (or  first)  w^  of  that  section  and  used  in  the  parabolic  fit. 

For  purposes  of  determining  the  parabola,  the  value  of  bound  vorticity  on 
the  extra  strip  is  taken  as  equal  to  the  value  at  the  midspan  of  the  last 
ordinary  lifting  strip  of  the  section,  even  though  this  is  not  strictly  true 
unless  the  slope  of  the  bound  vorticity  on  the  l3st  ordinary  strip  is  zero. 

7.12  The  Linear  Equations  for  the  Values  of  Surface  Source  Density 

A  dot  product  is  taken  of  each  source  velocity  V. •  at  each  on-body 

*  3 

control  point,  i  =  1,  2,  ....  N,  with  the  unit  normal  vector  of  the  surface 
element  containing  the  control  point.  Specifically, 

A..  =  n  •  ?  1  =  **  2’  N  (7.12.1) 

1J  1  1J  j  =  1,  2,  ...,  N 

The  scalar  N  x  N  matrix  A^  represents  the  normal  velocities  at  the 
control  points  due  to  unit  values  of  source  density  0  on  the  elements. 
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The  quantities 

N 

Y  A.  -a.  i  =  1,  2,  ....  N  (7.12.2) 

L _ »  •  J  J 

j=i 

where  the  source  densities  c ^  are  as  yet  unknown,  are  the  normal  velocities 
at  the  control  points  due  to  the  entire  surface  source  distribution.  For  the 
usual  condition  of  zero  normal  velocity  at  the  control  points  (7.12.2)  must 
be  set  equal  to  the  negative  of  the  normal  velocities  due  to  the  onset  flow. 
This  is  done  for  each  onset  flow.  Normal  components  of  the  basic 
circulatory  onset  flows  (7.11.3)  or  (7.11.7)  are  obtained  by  taking  dot 
products  with  the  unit  normal  vectors  in  a  manner  similar  to  (7.12.1),  i.e., 

n\k)  =  n.  •  vjk)  k  =  1,  2,  ....  L  (7.12.3) 

where  L  is  the  number  of  lifting  strips.  The  same  is  done  for  the  uniform 
onset  flow  V  ,  i.e., 

oo 

nH  =  n. .  *  ?  (7.12.4) 

1  1  03 

As  discussed  in  section  6.7  more  than  one  uniform  onset  flow  may  be  considered 
simultaneously,  In  which  case  there  is  an  n)  '  for  each  of  them. 

The  linear  equations  that  yield  the  values  of  source  density  on  the 
elements  are 

N  i  =1,2 . N 

£  A..c<k)  =  -n\k)  k  =  1,  2,  ....  L,  -  (7.12.5) 

These  are  solved  by  a  direct  elimination  procedure.  There  is  a  set  of  N 
values  of  a  ■  for  each  onset  flow,  including  all  uniform  onset  flows. 

7.13  Application  of  the  Kutta  Condition 

For  each  uniform  onset  flow  a  single  combined  set  of  source  densities 
is  calculated  from 

L 

Oj  =  c\m)  +  y  B(k)c^k)  j  =  1,  2,  ...,  N  (7.13.1) 

k=l 
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( k) 

where  L  is  the  number  of  liftinci  strips  end  where  the  B  1  are  as  yet 

(k) 

unknown.  The  combination  constants  B'  '  are  the  values  of  the  streamwise 
dipole  derivative  (bound  vorticity)  on  the  lifting  strips.  Similarly  there 
is  a  single  combined  onset  flow 


V^0)  =  ^(c)  +  y  B(k)^(k) 
1  1  L _ j  1 

k=l 

The  total  velocity  at  any  point  is 

fi  ‘  Z  Vo  '■ 

j=l 

=  +  £  B(k^k) 

k=l 

where  the  velocities 

N 

_  V  w-  (°°)  +  rr 

■  L->  ij  j 

j=l 

N 


V  ■  Z  V^’  +  ^ 


i  5  1,  2,  N  +  0 


i  =1,2,  . . . ,  N  +  0 


k  =  1 ,  2,  . . . ,  L 


(7.13.2) 


(7.13.3) 


(7.13.4) 


are  the  velocities  at  the  control  points  for  the  individual  onset  flows.  It 
is  important  to  point  out  that  velocities  (7.13.4)  are  calculated  only  for  the 
points  where  the  Kutta  condition  is  to  be  applied.  Only  the  velocity  (7.13.3) 
is  evaluated  at  all  points. 


As  mentioned  in  section  6.5,  there  are  two  rather  different  means  of 
applying  the  Kutta  condition. 

7.13.1  Flow  Tanqency  in  the  Wake 

The  first  method  for  applying  the  Kutta  condition  is  based  on  property 
(a)  of  section  6.5,  i.e,,  the  condition  that  a  stream  surface  of  the  flow 
leave  the  body  at  the  trailing  edge.  This  is  implemented  by  inputting  L 
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points  and  L  normal  vectors.  The  prints  are  considered  to  be  the  first  L 
off-body  points,  and  both  points  and  normal  vectors  are  designated  by  subscripts 
i  =  N  +  1 ,  N  +  2,  . . . ,  N  +  L .  The  total  velocity  at  these  points  is  given  by 
U.ld.d)  tor  these  values  of  i.  The  dot  product  of  each  velocity  is  taken 
with  the  corresponding  input  normal  vector,  which  is  presumed  to  be  the  unit 
normal  vector  tc  the  stream  surface.  The  results  are  set  equal  to  zero,  i.e., 


L 


k=l 


N  +  1, 


N  +  L  (7.13.5) 


Thus,  there  are  L  linear  equations  for  the  L  unknown  values 


,  namely 


L 


E 


DikB 


00 


-D. 

1 


where 


[)..  -  n.  •  ^ 

i  k  l  i 

D  •  =  n.  •  r.  ' 

i®  i  i 


i  =  N  +  1 ,  . . . ,  N  +  L  (7.13.6) 


k  =  1,  2, 
i  =  N  +  1 


L 

(7.13.7) 

,  N  +  L 


If  more  than  one  uniform  onset  flow  is  considered,  the  same  matrix  applies 

to  all  of  them.  Only  the  D.  are  different. 

|  oo 


7,13.2  Pressure  Equality  on  Upper  and  Lower  Surface  at  the  Trailing  Edge 


The  second  method  of  apolying  the  Kutta  condition  is  based  on  property 
(b)  of  section  6.5,  i.e.,  the  condition  that  the  pressures  be  equal  at  the  two 
control  points  of  each  strip  that  are  adjacent  to  the  trailing  edge.  The 
pressure  at  any  point  is  uniquely  determined  by  the  square  of  the  velocity 
magnitude,  which  is 


I2.  -  ^  •  V.  -  *  v^))  +  2  £  (7f }  •  ?|k))B(k) 

k=l 

+  y^(v|k)  •  ^.m))B(k)B{ni)  =  +  2  Y 


k^l  m-1 


L  I. 


E I  "ik™6"''1 


(k)K(m) 


(7.13.8) 


k=l  m=l 
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where  the  M's  are  defined  by  equation  (7.13.8).  Now  let  the  integer  q 
denote  the  lifting  strip,  i.e.,  q  =  1,  2,  ...,  L,  and  define 


!!  km  =  M ^ km  cont’ro^  Point  adjacent  to  trailing  edge  of  g-th  strip 
on  upper  surface)  ^  ^ 

-Mik,n  (at  control  point  adjacent  to  trailing  edge  of  q-th  strip 
on  lower  surface) 

Simi larly  define 


Vk  =  Mi»k  (upper  q-th)  “Mi-k  0ower  q‘th) 

H  =  (upper  q-th)  —  (lower  q-th) 


(7.13.10) 


where  tne  expressions  in  parentheses  in  (7.13.10)  are  intended  to  be 
abbreviations  of  the  Parentheses  in  (7.13.9).  With  this  notation,  the  equal- 
pressure  condition  is 


L  L 


"rZI  V«.B<k,B<"')  * 2  I  VkB<k)  +  V  =  0 


k=l  m=l 


q  =  1,  2,  ...,  L 


(k) 

This  represents  L  quadratic  equations  in  the  L  unknown  values  of  Er  . 
The  method  of  solution  is  a  Newton-Raphson  iterative  procedure.  Define  the 
derivative 


q  =  1,  2,  ....  L  (7.13.12) 
k  -  1 ,  2,  . . . ,  L 


Then  (7.13.11)  Is  solved  iteratively  by  solving  successive  sets  of  linear 
equations  for  the  changes  in  the  values  of  B^.  Namely, 


L 

GqkAB(k)  =  -Pq  a  -  1,  2,  ....  L  (7.13.13) 

k=l 

(k) 

At  any  stage  G„.  and  P„  are  evaluated  using  the  Bv  '  from  the  previous 

qk  q  (k)  (k) 

iteration.  Then  (7.13.13)  is  solved  and  new  B  '  computed  by  adding  aB  ' 

to  the  previous  values.  The  rate  of  convergence  of  this  process  nr  even  the 

existence  of  convergence,  cannot  be  proven  on  theoretical  grounds.  However, 
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in  virtually  all  cases  convergence  of  this  Iterative  process  has  been  very 
rapid.  There  can  be  difficulties,  however,  in  extreme  cases  (see  section  8.8). 
If  difficulties  should  arise  in  the  future  perhaps  (7.13.11)  should  be  solved 
by  a  different  iterative  procedure  than  that  represented  by  (7.13.13).  In  any 
event  the  procedure  of  section  7.13.1  can  be  used  with  confidence  since  no 
iteration  is  involved. 

If  several  uniform  onset  flows  are  considered,  the  same  H  ^  applies 
to  all  of  them. 


7.14  Final  Flow  Computation 


( k) 

Once  the  B  '  are  known,  a  single  set  of  source  densities  (for  each 
uniform  onset  flow)  is  computed  from  (7.13.1)  and  a  single  onset  flow  from 
(7.13.2).  Then  flow  velocities  at  the  on-body  control  points  and  off-body 
points  are  computed  from  (7.13.3).  Pressure  coefficients  at  control  points 
are  computed  from 


Forces  and  moments  are  computed  by  assuming  the  pressure  to  be  constant  ovc-r 
each  element.  If  denotes  the  area  of  the  i-th  element,  the  force  on 
this  element  Is 


^1  =  _niCp1Si 

and  the  moment  of  the  force  on  the  element  is 


m  *  "i  *  r( 


(7.14.2) 


(7.14.3) 


where  r^  represents  the  vector  displacement  of  the  control  point  of  the  ele¬ 
ment  from  an  Input  moment  reference  point.  With  the  above  assumption  forces 
and  moments  on  the  body  are  obtained  by  simple  summation 


(7.14.4) 
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Various  ranges  of  summation  are  used  in  (7.14.4)  so  that  forces  and  moments 
on  different  parts  of  the  configuration  may  be  calculated.  In  particular 
(7.14.4)  Is  performed  for:  each  nonlifting  section;  each  lifting  strip;  each 
lifting  section;  and  all  elements  of  the  entire  case. 

7.15  Computation  Time,  Effort,  and  Cost 

In  the  past  when  computing  machines  executed  one  program  at  a  time, 
computation  time,  effort,  and  cost  had  definite  and  agreed-upon  meanings.  The 
total  elapsed  time  necessary  to  execute  the  program  was  measured,  and  this  was 
charged  to  the  user  at  a  rate  of  a  certain  amount  of  money  per  hour.  Thus, 
computation  time  and  cost  were  simply  proportional.  Computational  effort  was 
slightly  less  direct,  since  elapsed  time  included  all  necessary  inputs  and 
outputs  and  certain  other  operations  in  addition  to  straightforward  arithmetic 
and  logic.  Nevertheless,  it  was  customary  simply  to  define  computational  effort 
as  the  time  required  to  execute.  Thus,  program  descriptions  customarily 
reported  computing  times,  but  by  implicit  assumption  they  were  also  defining 
computational  effort  and  cost. 

The  situation  was  changed  considerably  by  the  widespread  use  of  computer 
systems  that  process  several  unrelated  programs  simultaneously.  Computing 
time,  effort,  and  cost  are  no  longer  essentially  Identical;  and  indeed  their 
precise  relationship  cannot  be  specified,  except  possibly  in  terms  of  a 
particular  computing  facility.  Generally,  the  time  the  so-called  central 
processing  unit  spends  on  a  particular  program  is  recorded.  This  "CPU  time" 
is  that  required  for  the  arithmetic  and  logic  of  the  program.  From  CPU  time 
an  imaginary  "computing  time"  is  calculated  by  an  arbitrary  formula  that 
accounts  for  the  number  of  inputs  and  outputs.  Finally,  cost  is  determined 
by  multiplying  "computing  time"  by  a  rate  that  depends  in  a  complicated  way 
upon  the  fraction  of  the  total  capacity  of  the  computer  that  is  engaged  on 
that  particular  problem,  i.c.,  how  much  high-speed  core  storage  is  required, 
how  many  low-speed  tape  or  disk  units  are  used,  etc.  The  relationship  between 
CPU  time  and  "computing  time"  varies  from  facility  to  facility,  as  does  the 
formula  for  computing  cost  from  "computing  time."  Thus,  no  general  statements 
can  be  made.  A  change  in  the  accounting  procedure  can  significantly  alter  the 
cost  of  a  computer  run.  A  program  that  is  optimized  for  one  accounting 
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procedure  may  perform  poorly  on  another.  Often  the  use  of  less  high-speed 
storage  will  result  in  increases  in  computing  time  and  effort  but  a  decrease 
in  cost. 


While  nothing  definite  can  be  said,  still  there  is  a  need  for  some  simple, 
commonly-accepted  measure  of  the  size  of  a  program.  It  has  become  fairly 
common  to  use  CPI)  time  for  this  measure.  There  are  many  valid  objections  to 
this,  but  no  other  quantity  Is  more  acceptable.  It  should  always  be  •emembered 
that  CPU  time  is  merely  a  rough  guide  to  the  order  of  magnitude  of  the  program. 
For  the  present  application  CPU  times  are  given  for  the  IBM  370-165  computer. 

Below  are  CPU  times  obtained  for  typical  cases,  all  of  which  tiad  one 
plane  of  symmetry,  which  was  accounted  for  in  the  calculations.  The  element 
number  fl  refers  to  those  describing  one  half  of  the  body. 


Element 
Number  N 

250 

500 

650 

950 


CPU  Time 
in  Minutes 

1.7 

6 

12 

30 


The  times  for  the  lower  element  numbers  are  quite  acceptable.  The  rapid 
increase  In  CPU  time  with  element  number  for  the  larger  cases  is  presumably 
due  to  the  use  of  a  direct  solution  for  solving  the  simultaneous  equations. 
Clearly  an  iterative  solution  should  be  used  for  N  >  1000,  and  probably  for 
N  I  800.  On  the  other  hand,  the  direct  solution  is  seen  to  be  very  efficient 
for  N  <  500  and  probably  should  be  used  for  N  <  700. 
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8,0  NUMERICAL  EXPERIMENTS  TO  ILLUSTRATE  VARIOUS 
ASPECTS  OF  THE  METHOD 


8.1  Element:  Number  on  an  Isolated  Lifting  Wing 

It  is  important  in  three-dimensional  problems  to  he  able  to  estimate 
element  numbers  that,  make  the  error  in  the  potential  flow  calculation  consistent 
with  the  errors  inherent  in  the  approximation  of  a  real  flow  by  a  potential  flow, 
e.g.,  errors  due  to  neglect  of  compressibility  or  viscosity.  Too  small  an 
element  number  may  give  useless  results,  while  too  large  an  element  number 
leads  to  a  needlessly  large  computing  time.  For  good  accuracy,  complicated 
three-dimensional  geometries  require  more  elements  than  any  program  makes 
available  and  would  require  very  long  computation  times.  For  such  cases  the 
decision  regarding  element  number  is  an  easy  one;  simply  use  the  maximum  per¬ 
missible  number  of  elements  and  accept  a  lesser  accuracy.  For  simpler  cases 
a  study  of  the  matter  may  prove  worthwhile.  In  the  course  of  developing  the 
present  method  some  studies  of  this  type  were  conducted.  The  results  are 
included  here  in  the  hope  that  they  will  be  of  value  to  future  users.  Obviously, 
only  a  few  coses  could  be  studied  in  detail.  If  a  design  application  involves 
many  cases  of  similar  geometry,  an  element-number  study  for  that  particular 
geometry  should  be  conducted  by  the  user  before  proceeding. 

The  simplest  case  Is  that  of  an  isolated  wing.  Two  questions  must  be 
answered.  Mow  many  lifting  strips  should  be  placed  across  che  span  of  the 
wing?  How  many  lifting  elements  should  lie  on  each  strip?  The  second  of 
these  questions  can  be  answered  by  running  two-dimensional  cases  using  tne 
method  of  reference  1.  These  cases  are,  of  course,  very  fast  and  cheap 
compared  to  the  three-dimensional  cases  that  must  be  run  to  answer  the  first 
question.  For  this  investigation,  as  well  as  some  others  to  be  discussed 
below,  the  geometry  chosen  was  an  untwisted  wing,  which  is  described  fully  in 
reference  12.  The  planfom.  is  shown  in  figure  25,  and  the  airfoil  shape  in 
sections  parallel  to  the  symmetry  plane  of  the  wing  is  symmetric  ar.d  is  7.C  < 

percent  V  w  x.  Two-dimensional  considerations  lead  to  the  use  of  30  lifting 
elements  on  each  strip  -  15  on  each  of  the  upper  and  lower  surfaces.  This 
appears  to  be  aoout  a  minimum  number  for  acceptable  accuracy,  but  on  the  other 
hand  it  appears  sufficient  for  most  applications. 
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Calculations  were  performed  for  this  wing  with  various  numbers  of  lifting 
strips.  Four  of  the  cases  are  shown  in  figure  25.  They  range  from  6  to  20 
lifting  strips  on  the  right  half  of  the  wing.  In  comparing  solutions  the 
quantity  used  Is  the  local  section  lift  coefficient  as  a  function  of  spanwlse 
location.  This  quantity  Is  obtained  from  a  numerical  Integration  of  the  cal¬ 
culated  pressure,  which  Is  assumed  to  be  constant  over  each  surface  element. 

As  explained  In  section  9.1,  this  quantity  Is  considerably  more  sensitive  than 
pressure  distribution  in  the  sense  that  two  pressure  distributions  that  appear 
nearly  identical  may  have  section  lift  coefficients  that  are  noticeably  dif¬ 
ferent,  but  the  reverse  is  never  true.  The  cases  run  to  investigate  the 
effect  on  accuracy  of  the  number  of  lifting  strips  used  the  “step  function" 
bound  vorticity  option  (section  6.3)  and  applied  the  Kutta  condition  by  means 
of  the  condition  of  equal -pressure  at  the  first  and  last  control  points  of 
each  lifting  strip  (section  7.13.2).  Calculated  section  lift  coefficients  at 
eight  degrees  angle  of  attack  are  shown  in  figure  26  for  cases  of  8,  13,  and 
20  strips  (figure  25).  The  results  for  13  and  20  strips  are  nearly  identical 
except  for  a  small  region  near  90  percent  semispan,  and  the  20-strip  results 
are  thus  taken  as  correct.  The  values  of  lift  calculated  for  8  strips  are 
somewhat  too  large  but  may  be  close  enough  for  many  purposes.  However,  it 
appears  that  If  13  strips  are  used,  accurate  results  are  obtained,  and  this 
Is  thus  the  recommended  neighborhood  for  the  number  of  lifting  strips.  Thus, 
In  the  present  example  a  total  of  30  times  13  or  390  lifting  elements  are 
desirable. 


8.2  Two  Forms  of  the  Kutta  Condition 

In  section  5.5  two  forms  of  the  Kutta  condition  are  described.  They  may 
be  denoted  toe  wake- tannery  condition  (property  (a)  of  section  6.5)  sr$a  the 
eowfl' -ore* sure  condition  (property  (b)  of  section  6.5).  In  figure  15  calcu¬ 
lated  results  are;  compared  for  a  two-dimensional  case  where  the  stream  surface 
vl the  body  is  known  lx.-  lie  along  the  traillng-edge  bisector.  For  a  w'ng 
?f  the  type  shown  in  figure  25,  the  theory  of  reference  11  (section  6.5  and 
figure  53)  states  the  stream  surface  leaves  the  wing  along  the  tangent 
to  the  ;p,‘>er  surfeca.  However,  as  discussed  in  section  6.5  and  reference  6, 
it  is  often  more  accurate  to  apply  the  wake-tangency  condition  along  the 
twill ir-c-edge  bisector.  The  8- strip  case  of  figure  25  was  run  at  8  degree 
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angle  of  attack  using  the  "step  function"  option  for  bound  vorticity  with  a 
wake-tangency  condition  applied  at  a  distance  of  2  percent  of  local  chord 
from  the  trailing  edge.  Calculated  section  lift  coefficients  are  shown  in 
figure  26  for  points  of  application  of  the  wake-tangency  condition  lying  on 
the  traillng-edge  bisector  and  also  on  the  upper-surface  tangent.  For  the 
8-strip  case  the  error  •‘or  the  case  where  the  traillng-edge  bisector  is  used 
Is  seen  to  be  about  twice  as  large  as  that  obtained  with  the  equal-pressure 
condition  out  to  about  80-percent  semispan.  Application  of  the  wake-tangency 
condition  at  a  point  on  the  tangenc  to  the  upper  surface  gives  results  that 
are  very  seriously  in  error. 

Based  on  the  above  results  the  equal  pressure  condition  appears  superior 
to  the  wake  tangency  condition,  for  ordinary  cases.  Unless  otherwise 
Indicated,  It  is  used  for  all  cases  presented  In  this  report. 

8.3  Step  Function  and  Piecewise  Linear  Bound  Vorticity 

As  discussed  in  section  6.3,  the  present  method  has  two  options  for 
treating  the  variation  of  the  bound  vorticity  over  the  small  but  finite  "span" 
of  a  lifting  strip.  The  bound  vorticity  may  be  taken  either  constant  or 
linearally  varying  over  the  "span"  of  each  strip  to  yield  an  overall  spanwise 
variation  over  the  wing  that  is,  respectively,  a  step  function  or  a  piecewise 
linear  function  (figure  10).  To  investigate  the  differences  between  these 
two  representations  of  the  bound  vorticity,  the  13-strip  wing  of  figure  25 
was  run  at  8  degrees  angle  of  attack  with  an  equal-pressure  Xutta  condition 
using  each  of  the  bound-vorticity  options.  For  e  Th  case  the  bound  vorticity 
as  a  function  of  "spanwise"  location  on  the  wing  wh  oDl-’ned  by  fairing  a 
smooth  curve  through  the  computed  values  of  bound  vorticity  at  the  "midspans" 
of  the  lifting  strips.  Thus,  in  comparing  the  bound  vorticity  functions 
computed  by  the  two  options,  the  detailed  variation  over  the  individual  strips 
is  ignored.  The  calculated  results  are  shown  in  figure  27.  (Because  of  the 
sign  convention  adopted,  bound  vorticity  leading  to  a  positive  lift  has  a 
negative  value  of  the  proportionality  constant  B,  if  the  N-llne  is  Input 
with  the  lower  surface  first  as  recommended  In  sections  7.31  and  8.4.)  The 
results  are  seen  to  be  virtually  identical.  Surprlslrgly,  agreement  is  best 
In  region  of  rapid  variation  near  the  tip  and  worst  in  the  region  of 
relatively  slow  variation  near  the  plane  of  symmetry  of  the  wing. 
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To  further  compare  the  two  bound-vorticlty  options,  section  lift  coef¬ 
ficients  were  computed  by  numerical  Integration  of  the  surface  pressures.  The 
results  are  shown  In  figure  28.  Agreement  of  the  two  calculations  Is  good 
except  for  the  region  near  the  tip.  A  comparison  with  the  presumably  more 
accurate  results  from  the  20-strip  case  (figure  26)  Indicates  that  the  section 
lift  coefficients  calculated  by  the  step  function  option  are  more  accurate 
than  those  calculated  by  the  piecewise  linear  option.  The  values  of  pressure 
near  the  tip  are  affected  by  the  spanwise  velocity  component,  which  Is  sensi¬ 
tive  to  the  details  of  the  parabolic  fit  used  at  the  wing  tip  to  extrapolate 
the  piecewise  linear  bound  vortlcity  to  a  zero  value  In  the  fluid  (sections 
6.3  and  7.11).  However,  a  limited  amount  of  experimentation  with  the  para¬ 
bolic  fit  failed  to  Improve  the  calculated  distribution  of  section  lift 
coefficient  near  the  tip. 

Based  on  the  above  results  it  Is  concluded  that  there  is  no  apparent 
advantage  to  using  the  more  complicated  piecewise  linear  form  of  the  bound 
vortlcity,  at  least  for  simple  cases.  Accordingly,  the  simpler  step  function 
form  of  the  bound  vortlcity  has  been  used  for  all  cases  presented  in  this 
report.  However,  further  experimentation  with  the  piecewise  linear  form  of 
the  bound  vortlcity  seems  to  be  desirable,  particularly  for  more  complicated 
geometries.  Evidently,  an  Improved  wing  tip  condition  would  be  desirable. 

8.4  Order  of  the  Input  Points 

As  discussed  in  section  7.3.1  the  Input  can  be  arranged  so  that  the 
points  on  an  N-line  are  Input  In  one  of  two  orders.  In  any  case  the  first 
point  Input  is  at  the  trailing  edge.  Then  the  points  may  be  input  along  the 
lower  surface  of  the  wing  to  the  leading  edge  and  back  along  the  upper  surface 
to  the  trailing  edge.  Alternatively,  the  points  may  be  Input  along  the 
upper  surface  to  the  leading  edge  and  back  to  the  trailing  edge  along  the 
lower  surface.  (Recall  that  a  different  order  for  the  N-llnes  Is  required 
in  each  case.)  The  distinction  between  these  two  cases  is  Illustrated  in 
figure  22.  It  Is  concluded  in  section  7.3.1  that  the  calculated  values  of 
bound  vortlcity  should  be  equal  (corresponding  proportionality  constants  B 
equal  In  magnitude  and  of  opposite  sign)  In  the  two  cases  and  that  the  differ¬ 
ence  between  the  two  calculated  results  (figure  22c)  should  vanish  In  the 
limit  of  Infinite  element  number. 
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The  situation  described  above  was  Investigated  by  calculating  flow  about 
the  8-strlp  wing  of  figure  25  at  8  degrees  angle  of  attack  using  both  possible 
orders  for  the  Input  points.  Both  cases  used  the  step  function  option  for 
bound  vorticlty  and  applied  the  Kutta  condition  by  means  of  the  equal- 
Dressure  condition.  Calculations  were  performed  using  an  "open"  wing  tip 
>.  finite  thickness  and  repeated  using  a  "closed"  wing  tip,  for  which  the 
section  curve  at  the  tip  was  arbitrarily  given  zero  thickness.  There  was  no 
essential  difference  between  results  for  the  open  and  closed  tips,  so  only 
the  former  case  Is  presented  here.  Figure  29  compares  calculated  spanwise 
bound  vorticity  distributions  obtained  for  the  two  orders  of  input  points. 

The  two  distributions  are  seen  to  be  virtually  identical,  as  predicted. 

Figure  30  compares  calculated  spanwise  distributions  of  section  lift  coef¬ 
ficient,  which  are  obtained  by  integrating  surface  pressures.  Agreement  is 
good  except  near  the  wing  tip  where  the  solution  obtained  by  inputting  the 
lower-  surface  first  is  clear 1y  to  be  preferred.  What  has  occurred  is  that 
the  difference  of  the  two  solutions,  represented  by  the  solution  of  figure  ?2c, 
does  not  vanish  near  the  tip  because  of  the  finite  element  number. 

On  the  other  hand,  effects  like  that  of  figure  30  do  not  always  occur. 

Two  of  the  wing-fuselages  of  sect's"  9.3  were  computed  using  an  order  of  input 
points  such  that  the  upper  surface  of  each  section  curve  of  the  wing  was  input 
before  the  lower  surface.  Moreover,  the  wing  tips  in  both  cases  were  of  the 
"open"  type.  The  calculated  spanwise  distributions  of  section  lift  coefficient 
(figures  40a  and  42a)  appear  reasonable.  Of  possible  importance  is  the  fact 
that  the  strip  of  elements  adjacent  to  the  wing  tip  Is  considerably  wider  in 
both  wing-fuselage  cases  than  in  the  8-strip  winy  of  figure  25.  Evidently 
this  matter  deserves  further  study.  However,  inputting  the  lower  surface 
first  has  never  led  to  trouble. 

It  Is  concluded  that  ordering  the  input  so  that  the  lower  surface  ;jf  a 
lifting  section  is  input  before  the  upper  surface  is  a  desirable  procedure, 
and  It  is  followed  In  all  cases  presented  in  this  report  unless  otherwise 
stated.  The  terms  "lower"  and  "upper"  refer  to  the  usual  case  of  a  wing  at 
positive  angle  of  attack.  The  essential  condition  is  the  orientation  of  the 
surface  to  the  direction  of  the  onset  flow.  Thus,  for  a  general  flow  the 
term  "lower"  should  be  replaced  by  "windward"  and  the  term  "upper"  by 
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"leeward."  If  In  any  application  there  is  difficulty  deciding  which  side  of 
a  lifting  body  is  leeward  and  which  windward,  then  almost  certainly  it  will 
make  little  difference  which  is  chosen.  Finally,  the  differences  in  the 
calculated  results  for  the  two  orders  of  input  are  small  except  near  a  wing 
tip. 

8.5  Location  of  the  Trailing  Vortex  Wake 

As  discussed  in  section  6.3,  the  location  of  the  trailing  vortex  wake 
must  be  furnished  as  input  to  the  program.  In  practical  applications  the 
exact  location  is  not  known,  but  an  approximation  may  be  estimated  from 
experience.  To  determine  the  sensitivity  of  the  calculated  results  to 
wake  location,  several  geometries  were  calculated  with  different  wake  loca¬ 
tions.  Among  the  geometries  considered  was  the  wing  described  in  section  8.1 
and  another  wing  of  Identical  planform  with  camber  and  twist.  Wakes  were 
assumed  that  left  the  trailing  edge  along  the  bisector  and  also  along  the 
tangent  to  the  upper  surface.  Straight  wakes  were  used  and  also  wakes  that 
curved  and  became  parallel  to  the  direction  of  the  uniform  onset  flow.  None 
of  these  permutations  gave  any  significant  change  in  the  surface  pressures 
or  lift  distributions  on  the  wing.  It  is  thus  concluded  that  for  ordinary 
moderate  values  of  angle  of  attack,  trailing  edge  angle,  and  degree  of  camber 
any  reasonable  wake  location  gives  a  satisfactory  solution. 

It  may  be  recalled  that  the  two  solutions  of  figure  26  obtained  for 
a  wake-tangency  type  of  Kutta  condition  differed  very  markedly  from  each  other. 
This  was  due  to  the  locations  of  the  point  of  application  of  the  wake-tangency 
condition  not  to  the  assumed  wake  location. 

As  part  of  the  present  study,  a  review  of  the  literature  on  wake  location 
was  carried  out.  In  view  of  the  above,  the  results  of  the  review  do  not 
appear  to  be  of  paramount  importance  to  the  present  method.  This  is  fortunate 
because  the  amount  of  published  information  on  this  subject  is  not  very  large. 
The  literature  review  is  summarized  In  Appendix  B. 

It  should  be  emphasized  that  what  was  proved  in  the  above  study  is  that 
the  flow  on  a  lifting  body  is  insensitive  to  the  position  cf  its  own  wake. 
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Obviously  If  the  wake  from  one  lifting  hody  passes  near  another  body,  the 
flow  on  the  second  body  is  sensitive  to  the  location  of  this  wake.  This 
occurs,  for  example,  In  problems  of  wing-tail  interference. 

8.6  A  Wing  in  a  Wall.  Fuselage  Effects 

A  very  common  application  of  the  present  method  is  a  wing-fuselage.  For 
an  isolated  fuselage  much  larger  surface  elements  can  be  used  to  obtain  good 
accuracy  than  can  be  used  for  a  wing.  The  question  then  arises  as  tc  whether 
this  same  rather  sparse  element  distribution  can  be  used  for  a  fuselage  on 
which  a  wing  is  mounted.  To  investigate  this  point,  calculations  were  per¬ 
formed  for  a  straight  wing  protruding  from  a  plane  wall.  The  basic  geometry 
is  shown  in  figure  31a.  The  wing  has  a  rectangular  planform  with  span  equal 
to  five  times  chord.  The  airfoil  section. is  a  symmetric  one  with  a  thickness 
of  10-percent  chord.  The  plane  wall  extends  a  distance  of  five  airfoil  chords 
from  the  airfoil  in  both  fore-and-aft  and  sideways  directions. 


Two  studies  were  performed.  In  both  of  them  the  uniform  onset  flew  is 
parallel  to  the  plane  of  the  wall  and  is  at  10-degrecs  angle  of  attack  with 
respect  to  the  wing.  In  the  first  study  the  width  of  the  "extra  strip"  of 
elements  that  lies  on  the  opposite  side  of  the  wall  from  the  wing  was  given 
a  fixed  span  equal  to  one  airfoil  chord  as  shown  in  figure  31a.  Three  ele¬ 
ment  distributions  on  the  wall  were  used,  as  shown  In  figures  31b,  31c,  and 
31d.  The  dense  element  distribution  of  figure  31b  has  wall  elements  of  the 
same  chordwise  extent  as  the  elements  on  the  wing,  while  the  sparse  distribu¬ 
tion  of  figure  31 d  has  only  two  wall  elements  over  the  span  of  the  wing. 
Section  lift  coefficients  on  the  wing  calculated  with  the  three  different  wall 
element,  distributions  differ  by  one  unit  in  the  fourth  decimal  place,  which 
Is  utterly  negligible. 

The  second  study  used  the  wall  element  distribution  shown  in  figure  31 d 
and  considered  three  different  spanwise  extents  for  the  "extra  strip:"  one 
chord,  as  shown  in  figure  31a,  three  chords,  and  one-third  chord.  Calculated 
values  of  section  lift  coefficients  on  the  wing  differ  by  one  unit  in  the 
second  decimal  place.  This  is  of  some  importance  but  a  very  large  range  of 
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spans  Is  being  considered.  Certainly  It  can  be  concluded  that  the  span  of 
the  extra  strip  is  not  crucial. 

The  spanwise  variation  of  section  lift  coefficient  at  10  degrees  angle 
of  attack  for  the  case  with  a  one-chord  extra  strip  is  compared  in  figure  32 
with  that  obtained  at  the  same  angle  of  attack  for  the  Isolated  wing  of  aspect 
ratio  5,  and  that  for  the  aspect  ratio  10  wing  obtained  by  reflecting  the 
wing  in  the  plane  of  the  wall.  This  last  case  corresponds  to  use  of  an 
infinite  plane  wall.  It  can  be  seen  that  the  wall  of  figure  31  has  almost 
the  same  effect  as  the  infinite  wall.  The  difference  lies  not  In  the  finite 
element  size  but  in  the  finite  extent  (5  chords)  of  the  wall. 

8.7  A  Sudden  Change  in  Element  Shape 

Section  9.3  presents  results  for  a  wing  of  rectangular  planform  mounted  a? 
a  mldwlng  on  a  rectangular  fuselage.  Section  10.1  investigates  the  effects 
of  external  stores  mounted  on  this  wing-fuselage.  As  part  of  this  latter 
study,  two  different  element  distributions  were  used  on  the  wing.  These 
distributions  are  shown  In  figure  33.  In  both  cases  the  spanwise  distribu¬ 
tion  of  lifting  strips  is  identical.  In  the  first  case  the  distribution  of 
elements  Is  identical  at  all  spanwise  locations  (input  point  distribution 
Identical  on  all  N-lines)  so  that  the  elements  are  all  rectangular  and  are 
distributed  "straight"  across  the  wing.  In  the  second  case  "slanted"  elements 
are  used  on  four  consecutive  strips  near  mldsemlspan  (point  distribution 
changed  on  three  consecutive  N-lines).  In  both  cases  all  input  points  are 
exactly  on  the  wing  surface.  Tiie  freestream  was  taken  at  6  degrees  angle  of 
attack.  Calculated  spanwise  variations  of  section  lift  coefficient  are  shown 
in  figure  34a.  The  sudden  change  in  element  shape  causes  a  noticable  "wiggle" 
in  the  calculated  spanwise  lift  distribution.  In  a  more  complicated  applica¬ 
tion,  such  an  effect  mignt  be  taken  as  physically  real.  Accordingly,  if 
element  distributions  must  change  over  a  body,  it  is  preferable  that  they 
do  so  gradually. 

Figure  34b  compares  calculated  chordwise  pressure  distributions  at  the 
midspan  of  one  of  the  two  central  strips  of  the  slanted-element  region.  It 
appears  that  differences  in  lift  are  due  almost  entirely  to  differences  in 
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pressure  In  the  neighborhood  of  the  upper-surface  pressure  peak.  Elsewhere 
the  two  calculated  pressure  distributions  agree  very  well. 

H .8  An  Extreme  Geometry 

The  numerical  experiments  of  the  previous  sections  provide  guidelines 
on  the  use  of  the  method  for  ordinary  design  applications.  To  delineate 
limits  cf  validity  of  the  method,  calculations  were  performed  for  a  case 
having  a  highly  deflected  flap  (figure  35).  As  may  be  inferred  from  the 
figure,  the  geometry  shown  is  a  partial-span  flap  on  a  complete  wing-fuselage 
configuration  (reference  13).  This  portion  of  the  configuration  contains  the 
essential  difficulty,  and  it  was  selected  for  study  rather  than  the  complete 
wing-fuselage  to  save  computing  time.  This  geometry  was  selected  as  an 
extreme  example.  Real  flow  about  such  a  body  is  not  even  approximately  a 
potential  flow.  In  the  tests  of  reference  13  the  flow  over  the  geometry  of 
figure  35  was  separated  even  if  area  suction  was  used  on  the  body  surface. 

When  calculations  were  performed  at  zero  angle  of  attack  with  the  equal- 
pressure  Kutta  condition,  the  iterative  procedure  of  section  7.13.2  diverged 
strongly.  This  is  the  only  case  to  date  where  this  failure  occurred*  The 
wake-tangency  Kutta  condition  of  section  7.13.1  was  applied  and  gave  a  reason¬ 
able  spanwlse  distribution  of  bound  vorticlty.  However,  the  pressures  at  the 
two  control  points  of  each  strip  adjacent  to  the  trailing  edge  were  not 
approximately  equal.  In  a  case  such  as  this  the  proper  location  for  the  trail¬ 
ing  vortex  wake  cannot  be  approximated  well  by  intuition.  Calculations  were 
performed  with  different  assumed  wake  locations,  and  significant  differences 
In  the  calculated  flow  were  obtained. 

Thus,  it  appears  that  the  present  method  can  calculate  flow  about  "normal" 
configurations  in  a  routine  fashion  but  that  there  are  limits  beyond  which 
some  care  is  required. 


*However,  see  section  10.4 
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9.0  COMPARISON  OF  CALCULATED  RESULTS  WITH  EXPERIMENTAL  DATA 


9.1  General  Remarks 

In  the  following  sections,  flow  quantities  calculated  by  the  present 
method  are  compared  with  experimental  data.  All  computations  follow  the 
recommendations  of  section  8.0.  In  particular,  the  step  function  option  for 
bound  vorticity  and  the  equal -pressure  Kutta  condition  are  used. 

Two  flow  quantities  are  compared:  the  section  lift  coefficient  as  a 
function  of  spanwlse  location  and  the  chordwise  pressure  distribution  at 
fixed  spanwise  location.  The  former  of  these  is  much  more  sensitive  than 
the  latter.  As  will  be  seen,  the  usual  situation  is  one  in  whicn  the  calcu¬ 
lated  and  experimental  pressure  distributions  agree  fairly  well  but  the  sec¬ 
tion  lift  coefficients  are  noticeably  different  because  the  difference  between 
the  two  pressure  distributions  is  of  constant  sign  and  its  Integrated  effect  is 
signi ficant. 

It  Is  well  known  that  for  unseparated  flow  the  effect  of  viscosity  is 
small  in  nonlifting  flow  but  is  quite  significant  in  flows  with  lift.  While 
the  exact  magnitude  of  the  effect  depends  on  the  Reynolds  number,  the  general 
effect  of  viscosity  Is  to  reduce  the  lift  about  10  percent  from  its  inviscid 
value.  In  two  dimensions  calculated  Inviscid  and  experimental  pressure  dis¬ 
tributions  on  an  airfoil  are  quite  different  If  they  correspond  to  equal 
angles  of  attack  but  agree  very  well  if  they  correspond  to  equal  lift  coef¬ 
ficients.  That  Is,  the  principal  effect  of  viscosity  is  on  the  lift  rather 
than  on  the  details  of  the  pressure  distribution.  This  last  is  probably 
true  In  three  dimensions  also.  However,  a  condition  of  equal  lift  is  difficult 
to  arrange  if  there  is  a  spanwlse  variation  of  the  lift  coefficient  and  of  the 
corresponding  viscous  effect.  In  any  case  it  is  desireable  to  calculate  the 
lift,  not  to  accept  it  as  given.  Thus,  the  proper  aim  is  to  calculate  cor¬ 
rect  flow  quantities  at  a  given  angle  of  attack.  Accordingly,  comparisons 
of  calculated  and  experimental  results  are  given  here  at  equal  angles  of 
attack.  It  is  believed  that  most,  if  not  all,  of  the  differences  between  the 
calculated  and  experimental  quantities  are  due  to  the  effects  of  viscosity 
(and  compressibility  in  some  of  the  tests).  Some  preliminary  work  on  this 
matter  has  been  done  and  confirms  this  opinion  (See  the  following  section). 
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9.2  An  Isolated  Wirtq 


An  untwisted  swept  wing  with  a  symmetric  airfoil  section  is  described 
in  section  8.1.  Low  speed  wind  tunnel  data  are  available  for  this  wing  in 
reference  12.  At  a  Reynolds  number  of  18  million  the  results  indicate  that 
no  separation  occurs  at  an  angle  of  attack  of  8  degrees,  and  calculations 
and  experiment  are  compared  for  this  flow  condition.  Results  are  shown  in 
figure  36.  It  can  be  seen  that  calculated  and  experimental  pressures  agree 
rather  well  at  all  chordwlse  and  spanwlse  locations,  except  possibly  near  the 
trailing  edge  near  the  tip  (figure  36d).  Calculated  and  experimental  distri¬ 
butions  of  section  lift  coefficient  are  quite  similar  in  shape,  but  the  cal¬ 
culated  inviscld  values  are  too  high  by  10-15  percent. 

To  test  the  hypothesis  that  viscous  effects  are  primarily  responsible 
for  the  disagreement  between  calculation  and  experiment,  a  crude  estimate 
of  the  distribution  of  boundary-layer  displacement  thickness  was  added  to 
the  wing.  Flow  about  the  altered  body  was  calculated  by  the  present  method, 
and  the  results  are  also  shown  in  figure  36.  A  dramatic  Improvement  in  the 
spanwlse  lift  distribution  is  evident  in  figure  36a.  Thus,  the  hypothesis 
concerning  viscous  effects  appears  valid.  Changes  in  the  pressure  distri¬ 
butions  are  less  spectacular,  but  as  mentioned  above,  these  are  relatively 
insensitive. 


9.3  Wing-Fuselages 

Reference  14  presents  experimental  data  for  a  simplified  wing-fuselage 
that  consists  of  an  uncambered  wing  of  rectangular  planform  mounted  as  a 
midwing  on  a  round  fuselage.  Low-speed  tests  were  conducted  at  the  very  low 
Reynolds  number  of  0.31  million.  Thus,  viscous  effects  are  rather  large  for 
this  experiment.  This  Is  not  a  very  suitable  case  for  comparison  with  a 
potential  flow  method.  It  was  selected  because  calculated  and  experimental 
results  for  a  very  similar  geometry  are  presented  in  reference  6,  and  it 
seemed  interesting  to  compare  the  predictions  of  the  present  method  with  that 
of  reference  6.  The  situation  is  complicated  by  the  fact  that  the  data  of 
reference  6  were  taken  at  a  higher  Reynolds  number  of  0.66  million,  so  that 
viscous  effects  are  reduced. 
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Figure  37  shows  the  geometry  of  the  configuration.  Figure  38  compares 
calculated  and  experimental  results  on  the  wing  for  an  angle  of  attack  of 
6  degrees.  The  two  spanwlse  lift  distributions  are  of  similar  shape  with  the 
calculated  values  about  20  percent  higher  than  the  experimental.  The  pressure 
distributions  are  in  better  agreement,  but  the  differences  in  lift  are  so 
great  that  the  pressures  on  the  upper  surface  are  affected.  No  conclusions 
can  be  drawn  regarding  the  relative  effectiveness  of  the  present  method  and 
that  of  reference  6.  The  agreement  of  calculation  and  experiment  presented 
In  reference  6  Is  much  the  same  as  that  shown  in  figure  38. 

A  configuration  of  current  interest  is  a  wing  with  a  so-called  "super¬ 
critical"  airfoil  section  mounted  as  a  high  wing  on  a  fuselage.  The  configu¬ 
ration  and  the  surface  elements  used  in  the  calculation  are  shown  in  figure 
39.  The  "supercritical"  airfoil  section,  which  is  also  shown  in  figure  39 
is  very  thin  in  the  neighborhood  of  the  trailing  edge  and  carries  a  relatively 
large  percentage  of  its  lift  in  this  region.  As  can  be  seen  in  figure  39, 
the  fuselage  represents  an  attempt  at  realism  with  low  element  numbers.  The 
cockpit  canopy  and  the  wing-tunnel  sting  are  both  accounted  for.  Figure  40 
compares  calculated  results  on  the  wing  at  7  degrees  angle  of  attack  with 
experimental  data  from  a  low-speed  wind-tunnel  test  conducted  by  Douglas 
personnel.  The  comparison  of  the  section  lift  coefficient  distributions 
exhibits  the  by-now-famillar  behavior  of  similar-shaped  curves  with  experi¬ 
mental  values  lower  than  calculated  ones  due  to  viscous  effects.  The  agree¬ 
ment  of  the  pressure  distributions  Is  quite  good,  especially  at  the  leading- 
edge  peak.  Also,  the  characteristic  "supercritical"  type  pressure  distribu¬ 
tion  aft  of  midchord  is  predicted  fairly  well  by  the  calculations.  The 
pressure  distributions  of  figure  40b  at  15  percent  semispan  are  at  a  location 
guite  near  the  wing-fuselage  junction,  which  Is  at  13.3  percent  semispan. 

Thus,  th.'ee-dimenslonal  interference  effects  are  relatively  large  at  this 
location  and  are  predicted  fairly  well. 

A  comparison  configuration  to  the  one  of  the  previous  paragraph  consists 
of  a  wing  with  a  conventional  airfoil  section  mounted  as  a  low  wing  on  a 
fuselage.  The  body  and  the  surface  elements  used  to  represent  it  are  shown 
in  figure  41.  Once  again  the  cockpit  canopy  and  the  wind  tunnel  sting  are 
accounted  for  in  trie  calculations.  Wind  tunnel  tests  of  this  configuration 
at  6.9  degrees  angle  of  attack  were  conducted  by  Douglas  personnel  at  a 
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freestream  Mach  number  of  0.5.  These  test  results  are  compared  with  the 
Incompressible  calculations  of  the  present  method  in  figure  42.  At  first  sight 
the  results  appear  quite  gratifying.  The  agreement  of  calculation  and  experi¬ 
ment  is  much  better  for  this  case  than  for  the  supercritical  wing-fuselage, 
whose  results  are  shown  In  figure  40.  Agreement  is  especially  good  for  the 
pressures  at  25-percent  semispan,  figure  42c,  but  the  span1  ise  distribution 
of  section  lift  coefficient  (figure  42a)  is  also  In  fairly  good  agreement. 
Unfortunately,  part  of  the  reason  for  this  agreement  is  that  the  errors  in  the 
calculation  due  to  neglect  of  viscosity  and  the  errors  due  to  neglect  of 
compressibility  are  of  opposite  sign  and  tend  to  cancel  each  other.  To 
illustrate  the  magnitude  of  the  compressibility  effect,  the  calculated  results 
have  beer  divided  by  the  quantity  Vl  —  M^,  where  M  denotes  freestream 
Mach  number  (figure  42).  This  type  of  correction  has  validity  in  two-dimensions 
within  the  limits  of  small  perturbation  theory,  hut  it  has  no  justification  in 
three-dimensions.  The  curves  with  this  divisor  in  figure  42  are  not  attempts 
to  quantitatively  predict  compressihility  effects  '  ct  are  supposed  to 
Illustrate  their  general  magnitude.  It  appears  that  when  compressibility  is 
accounted  for  the  agreement  of  calculation  and  experiment  for  the  configuration 
of  figure  41  Is  about  the  same  as  for  the  configuration  of  figure  39. 

The  discussion  of  the  previous  paragraph  also  points  out  the  need  for 
a  compressibility  correction  to  be  added  to  the  present  method.  Based  on 
previous  two-dimensional  experience,  this  should  prove  to  be  much  easier  than 
accounting  for  viscous  effects.  The  classical  procedut e  is  based  on  the 
Goethert  transformation.  However,  this  is  not  very  satisfactory.  Its  accuracy 
Is  poor  in  regions  such  as  wing  leading  edges  where  the  surface  slope  is  not 
approximately  parallel  to  freestream  velocity.  Moreover,  a  complete  calcula¬ 
tion  must  be  performed  from  the  beginning  for  each  Mach  number.  What  is 
needed  is  a  procedure  that  obtains  compressible  results  directly  from  an 
Incompressible  solution,  so  that  only  one  lengthy  flow  calculation  need  be 
performed  by  the  present  method.  An  example  of  such  a  method  is  presented 
In  reference  6,  but  the  results  are  not  entirely  satisfactory.  Evidently 
further  investigation  is  required. 

The  wings  of  the  configurations  of  figures  39  and  41  were  input  with  the 
upper  surface  first.  The  calculated  distributions  of  section  lift  coefficient 
that  are  given  in  figures  40  and  42  do  not  show  unusual  behavior  near  the  wing 
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lip, as  was  exhibited  In  figure  30.  These  are  the  cases  referred  to  in 
section  8.4. 

9.4  A  Wing-Fuselage  in  a  Wind  Tunnel 

A  rather  extensive  study  was  performed  for  the  somewhat  unusual  configura¬ 
tion  shown  In  figure  43.  Agreement  of  calculation  and  experiment  was  never 
obtained,  but  the  results  are  a  good  illustration  of  the  versatility  of  the 
method  and  the  uncertainty  connected  with  much  wind  tunnel  data.  The  basic 
configuration  Is  a  W-wIng  mounted  on  a  round  fuselage.  In  the  wind  tunnel 
the  model  was  mounted  on  a  support  strut,  as  shown  in  figures  43a  and  43b. 

Tne  data  <ere  supposedly  corrected  for  all  tunnel  interference  effects  (refer¬ 
ence  15).  Thus  the  initial  calculation  was  for  the  isolated  wing-body  (no 
strut  or  tunnel  walls)  at  a  corrected  free-air  angle  of  attack  of  4.4°,  which 
supposedly  corresponds  to  a  tunnel  angle  of  attack  of  4°.  A  comparison  of 
calculated  and  experimental  section  lift  coefficients  across  the  span  are 
shown  in  figure  44a.  Agreement  is  good  except  at  the  kink  and  near  the  tip, 
where  viscous  effects  are  important.  The  lack  of  response  of  the  calculations 
to  the  kink  was  surprising.  However,  an  approximate  potential  flow  calculation 
gives  results  that  agree  in  general  character,  but  not  in  precise  value,  with 
the  calculations  of  the  present  method.  This  also  Indicates  that  viscosity 
is  responsible  for  the  dip  at  the  kink  in  the  experimental  curve  of  lift  coef¬ 
ficient.  However,  the  agreement  of  calculated  and  experimental  pressure 
distributions  is  not  good,  as  is  shown  in  figures  44b  and  44c  for  two  spanwise 
locations,  A  check  on  the  blockage  and  upwash  corrections  that  were  applied 
to  the  data  raised  some  questions.  Accordingly,  calculations  were  performed 
for  the  strut-mounted  body  in  the  tunnel,  which  is  shown  In  figure  43.  The 
actual  wind  tunnel  angle  of  attack  of  4°  was  used.  The  results  of  this 
calculation  are  Included  in  figure  44.  Figures  44b  and  44c  show  a  rather  large 
effect  of  the  strut  on  the  lower  surface  pressures,  particularly  at  the 
inboard  location.  The  strut  effect  is  mainly  to  increase  blockage  below  the 
wing  but  not  above.  Thus,  lower-surface  pressures  are  lowered  (higher 
velocity)  and  the  result  is  the  loss  of  lift  shown  in  figure  44a.  The  effect 
shown  Is  much  larger  than  the  nominal  upwash  and  blockage  corrections,  and 
thus  some  doubt  exists  as  to  the  validity  of  the  data.  This  is  an  interesting 


application  of  the  program.  No  conventional  correction  could  account  for  the 
strut  effect,  but  the  program  obtains  it  rather  easily.  However,  the  experi¬ 
mental  pressures  are  still  more  negative  than  the  calculated  in  a  way  that, 
cannot  be  explained  on  physical  grounds.  A  difference  in  reference  static 
pressure  possibly  could  cause  this  discrepancy. 
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in. 0  INTERFERE?*.'*  STUDIES 


Presented  below  are  three  jxaijvples  '-f  the  use  o F  the  present  method  to 
predict  the  three* dimensionc!  lnterefereuce  effects  on  Sifting  wings  of  other 
todies  In  close  proximity  *  both  ‘‘Ifti^g  end  nonMfting.  The  cases  discussed 
‘p  section  1C.1  are  generated  c*or,rtriec  that  heve  all  the  essential  properties 
of  c-stual  designs,  However,  ic.<*y  are  not,  themselves,  deslons  of  Interest*  but 
se*v=-  tc  111  us  trite  the  capability  of  tne  nreseot  r'etned.  On  the  other  hand, 
the  cases  of  section;  10.2,  10,2,  and  10. *  were  generate*:  by  outside  users  and 
tn»?  results  thp-nscivec  were  of  Interest.  These  casus  thus  represent  the  first 
use  of  the  method  as  >.n  analysis  tool. 

* 

10.1  K1ng«Fuselt<5*-  with  External  -*tcr*s 

Fc:  this  exeraole  of  an  inter-torenc*:  study,  the  bas*c  geometry  *s  a 
rectangular  suing  nrounted  as  a  iHdvIng  on  a  round  ft-stlige.  This  geometry  Is 
shown  in  figure  V  anj  Is  discussed  in  section  9.3  hr*.?  reference  14.  Two 
externa ’-store  configurations  a»-e  considered.  The  first  consists  of  a  tip 
tank.  Thf  geometry  and  element  distributing  for  this  case  are  shown  In 
figure  Via.  Thy  secovJ  configuration  consists  of  the  same  external  store 
mounted  beneath  the-  wing  on  *  short  pylon  centered  at  6C-percent  senlspan. 

The*  seoretry  and  element,  distribution  for  this  case  are  shown  fn  figure  45b. 

rigurh  45fj  alsj  shows  the  "extra  strip"  of  elements  inside  the  tip  tank. 

?%  c'sc «;>sec  in  ioctton  3.8,  there  will  he  a  "hub  vortex"'  trailing  downstream 
f rcri!  the  tip  tank.  However,  this  does  net  appear  to  cause  any  numerical 
prob'e*?,  ano  calculated  surface  veioclties  near  the  downstream  end  of  the 
tip  tank  seev  ent^ely  reasonable,  Figure  450  snovrs  that  the  trailing  edge 
of  the  wing  is  continuous  across  the  span.  Accordingly,  the  "Ignored  element" 
oroc*dur*»  of  section  t.S  is  used  for  the  elements  on  the  lower  surface  of  the 
wing  that,  are  covered  or  partially  covered  by  the  pyloo,  as  illustrated  by 
the  dotted  planfonm  of  the  pylon  In  figure  45b. 

Calculations  were  performed  at  6  degrees  angle  of  attack  for  the  three 
configurations::  the  clean  wing,  the  wing  w*th  tto  tart,  and  the  wing  with 
oylon-sxyjnted  external  st*-" .  The  round  fuselage  was  present  In  all  cases. 


Calculated  results  on  the  wing  are  compared  in  figure  46.  Figure  46a  compares 
spanwlse  distributions  of  section  lift  coefficient.  The  addition  of  the  tip 
tank  to  the  wing  prevents  the  lift  from  falling  to  zero  at  the  tip,  so  there 
is  a  large  Increase  in  lift  coefficient  in  this  region.  Moreover,  as  pre¬ 
dicted  by  various  theories,  the  addition  of  a  tip  tank  Increases  the  effective 
aspect  ratio  of  the  wing  and  thus  increases  the  section  lift  coefficient  all 
the  way  to  the  fuselage.  The  effect  of  the  pylon-mounted  external  store  falis 
to  zero  at  the  wing  tip,  but  at  the  fuselage  this  effect  is  about  the  same 
size  as  that  of  the  tip  tank  but  in  the  opposite  direction.  The  major  effect 
of  the  pylon  is  to  reduce  lift  in  its  vicinity  by  increasing  lower-surface 
velocities  and  thus  reducing  lower-surface  pressures  (figure  46c).  Notice 
that  lift  on  the  wing  cannot  be  meaningfully  computed  at  the  spanwise  loca¬ 
tion  of  the  pylon  because  the  lower  surface  is  not  exposed  to  the  flow.  Of 
course,  there  Is  a  force  on  the  external  store,  but  it  cannot  be  meaningfully 
associated  with  a  particular  location  on  the  wing.  The  bound  vorticity 
distribution  is  continuous  across  the  span.  The  general  form  of  this  function 
Is  quite  similar  to  that  of  the  section  lift  coefficient.  Indeed,  it  looks  as 
If  a  human  had  faired  a  plausible  joining  curve  between  the  disjoint  portions 
of  the  curve  of  figure  46a. 

Chordwlse  pressure  distributions  for  the  clean  wing  and  for  che  wing  with 
the  tank  are  compared  in  figure  46b  for  a  spanwise  location  close  to  the  tip 
tank.  The  increase  in  lift  due  to  the  tip  tank  is  seen  to  be  primarily  due  to 
increased  velocity  on  the  upper  surface  of  the  wing.  Figures  46c  ana  46d 
compare  chordwise  pressure  distributions  for  the  clean  wing  and  for  the  wing 
with  pylon-mounted  external  store.  The  spanwise  location  of  figure  46c 
represents  the  strip  of  elements  ironediately  adjacent  to  the  pylon  location. 

The  considerable  reduction  in  lower  surface  pressures  due  to  the  presence  of 
the  pylon  and  store  is  evident.  Upper  surface  pressures  are  scarcely  affected. 
Figure  46d  compares  pressure  distributions  on  the  upper  surface  of  the  wing 
corresponding  to  a  strip  of  elements  that  lies  directly  above  the  location 
of  the  pylon.  The  pressure  distribution  computed  for  the  wing  with  pylon- 
mounted  external  store  is  quite  reasonable  except  for  a  "hump"  between 
66-percent  chord  and  90-percent  chord.  Examination  of  the  side  view  of 
figure  45b  shows  that  In  this  region  the  surface  elements  on  the  pylon  and  wing 
have  dimensions  that  are  considerably  larger  than  the  local  thickness  of  the 
wing.  Thus,  the  presence  of  pylon  is  sensed  "through"  the  wing  on  the  upper 
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surface.  An  Increase  in  element  number  could  remove  this  pressure  "hump"  but 
this  seems  unnecessary.  The  proper  way  to  fair  the  upper  surface  pressure 
distribution  is  quite  obvious.  It  is  felt  that  the  computed  results  of 
figure  46  represent  a  very  successful  application  of  the  present  method. 

10.2  Wing  with  Endplates 

A  case  in  which  the  calculated  results  were  of  interest  to  a  user 
concerned  the  effect  of  endplates  on  a  wing.  The  wing  in  question  has  a 
rectangular  planform  of  aspect  ratio  1.4  and  an  NACA  4415  airfoil  section. 

The  endplate  has  a  planform  consisting  of  a  semicircular  forward  section  and 
a  rectangular  rear  section.  The  entire  configuration  is  shown  in  figure  47. 
Three-dimensional  calculations  were  performed  at  10  degrees  angle  of  attack 
with  and  without  the  endolates.  A  two-dimensional  calculation  was  also 
obtained  for  comparison.  This  last  corresponds  to  a  case  of  endplates  of 
Infinite  extent. 

Calculated  "-esults  cro  compared  in  figure  48.  It  can  be  seen  from 
figure  48a  that  the  addition  of  the  endplates  produces  a  lift  distribution 
that  is  virtually  independent  of  spanwise  location.  (The  slight  drop  at  the 
last  spanwise  location  is  probably  a  numerical  error  and  should  be  faired  out. 
However,  the  level  of  the  lift  is  much  closer  to  that  of  the  isolated  three- 
dimensional  wing  than  it  is  to  the  two-dimensional  value.  The  chordwise 
pressure  distributions  in  the  symmetry  plane  (figure  48b)  also  exhibit  this 
behavior. 

In  performing  the  above  calculations  the  endplates  were  taken  as  simple 
symmetric  airfoils  4-percent  thick  and,  of  course,  had  sharp  trailing  edges. 

If  an  endplate  were  present  without  the  wing,  it  would  he  nonlifting..  In  the 
presence  of  the  wing  the  endplate  has  ar.  inward  lift  above  the  wing  and  an 
outward  lift  below  it.  The  level  of  lift  on  the  endplate  above  fhe  wing  is 
considerably  larger  than  that  on  the  endplate  below  the  wing  (about  three 
times)  and  is  about  one-fourth  the  level  of  lift  on  the  wing. 

This  case  differs  from  previous  cases  in  that  it  represents  an  inter¬ 
section  of  two  lifting  portions  of  a  configuration.  The  "ignore"  option  of 
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section  6,8  was  used  on  certain  strips  of  the  endplate  to  accommodate  the 
wing  Intersection, 

10,3  Wing  In  a  Wind  Tunnel 

The  wing  of  aspect  ratio  1.4  described  in  section  10.2  was  considered 
to  be  in  a  wind  tunnel  at  10  degrees  angle  of  attack,  as  shown  in  figure  49. 

If  the  wing  completely  spans  the  tunnel,  the  theoretical  inviscid  result  is 
the  two-dimensional  flow  about  the  airfoil  section  in  the  presence  of  the 
upper  and  lower  walls,  i.e.,  about  the  sldeview  of  figure  49  considered  as 
a  two-dimensional  flow.  However,  the  presence  of  the  gaps  between  the  wing 
tips  and  the  tunnel  sidewalls  allows  the  bound  vorticity  on  the  wing  to  fall 
to  zero  at  the  tips  and  introduces  significant  three-dimensional  effects. 

The  purpose  of  the  calculation  was  to  evaluate  these  three-dimensional  effects. 

Figure  50  compares  calculated  results  for  the  above-described  two- 
dimensional  case  with  those  for  the  three-dimensional  wing  with  and  without 
the  wind  tunnel  sidewalls.  All  cases  include  the  effects  of  the  top  and 
bottom  walls  of  the  wind  tunnel.  In  the  three-dimensional  case  without  side- 
walls,  the  top  and  bottom  walls  have  been  extended  horizontally  a  distance  of 
several  wing  spans.  The  importance  of  the  gaps  is  quite  evident  in  figure  50. 
Results  for  the  case  of  the  small  but  finite  gaps  are  much  closer  to  those 
for  infinite  gaps  (sidewalls  removed)  than  to  those  for  zero  gaps  (two- 
dimensional  case). 

10.4  Wing  With  Endp'iates  in  a  Wind  Tunnel 

As  a  final  example,  the  wing  with  endplates  (figure  47)  was  inserted 
in  the  wind  tunnel  shown  in  figure  49  to  obtain  the  configuration  shown  in 
figure  51.  When  calculations  were  performed  for  this  case  with  the 
equal-pressure  Kutta  condition,  the  iterative  procedure  of  section  7.13.2 
appeared  to  he  neutrally  convergent  and  thp  iterations  never  fully  "settled 
down".  This  may  have  been  caused  by  the  close  proximity  of  the  elements  on 
the  wind  tunnel  wail  to  the  trailing  edges  of  the  endplates.  In  all  cases 
except  this  one  and  the  strongly  divergent  case  of  section  8.8  the  iterative 
procedure  of  section  7.13.2  converged  very  rapidly. 


120 


Because  of  the  above  situation,  calculations  were  performed  for  the 
configuration  of  figure  51  using  the  wake-tangency  Kutta  condition.  Figure 
52  compares  the  results  obtained  with  those  for  the  isolated  wing  in  free 
air.  To  evaluate  the  effect  on  the  results  of  the  form  of  the  Kutta  condition, 
calculations  were  performed  for  the  wing  in  free  air  using  both  forms  of  the 
Kutta  condition.  As  can  be  seen  in  figure  52,  the  effect  of  the  form  of  the 
Kutta  condition  is  not  large,  and  most  of  the  differences  between  the  calcula¬ 
tions  for  the  wing  with  endplates  in  the  tunnel  and  the  various  other  results 
shown  in  figures  48,  50,  and  52  are  due  to  differences  between  the  geometries. 
It  is  evident  from  figure  52a  that  the  effects  of  endplates  and  wind  tunnel 
walls  together  give  a  lift  distribution  independent  of  spanwise  location. 
(Again,  the  drop  in  lift  at  the  spanwise  location  adjacent  to  the  endplate 
is  probably  a  numerical  inaccuracy  and  should  be  faired  out).  Moreover,  the 
level  of  the  lift  is  much  closer  to  the  two-dimensional  value  than  were  those 
obtained  using  endplates  or  tunnel  walls  separately  (figures  48a  and  50a). 

The  chordwise  pressure  distributions  of  figure  52b  also  show  the  inter¬ 
ference  effects  described  above.  Also  shown  are  the  small  but  noticeable 
differences  between  upper  and  lower-surface  trai ling-edge  pressures  in  the 
cases  that  used  the  wake-tangency  Kutta  condition.  For  the  wing  in  free 
air  the  pressure  distributions  calculated  using  the  two  forms  of  the  Kutta 
condition  differ  from  each  other  only  'in  the  vicinity  of  the  trailing  edge 
and  are  essentially  identical  over  most  of  the  surface. 
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APPENDIX  A 

RELATION  BETWEEN  DIPOLE  ANO  VORTEX  SHEETS  OF  VARIABLE  STRENGTH 


Consider  a  surface  S  in  space  bounded  by  a  closed  curve  c.  (If  S 
is  a  closed  surface,  c  vanishes.)  At  any  point  U,  n>  O  on  S  the 
unit  normal  vector  is  n,  and  at  any  point  on  c  the  unit  tangent  vector 
is  t.  The  vector  between  the  point  (5,  n.  c)  and  a  general  point 
(x,  y,  z)  In  space  Is  denoted  if,  and  the  length  of  this  vector  is 
denoted  R.  Specifically, 


R  =  (x  -  c)T  +  (y  -  n)3*  +  (z  - 
R  =  V(x  -  c)2  +  (y  -  r,)2  +  (z  -  o2 


(A-l) 


The  gradient  operator  grad^  is  used  to  denote  that  derivatives  are  taken 
with  respect  to  x,  y,  z.  Similarly,  grad  is  the  gradient  operator  that 
differentiates  with  respect  to  Ci  n»  c. 


THEOREM:  Let  the  surface  S  be  covered  with  a  variable  dipole  distri¬ 
bution  of  intensity  v .  (The  dipole  axes  are  along  n)  The  velocity  at 
(x,y, z)  due  to  the  dipole  sheet  is  equal  to  the  sum  of  the  velocities  due  to 
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a  certain  vortex  sheet  of  strength  to  on  S  and  due  to  a  vortex  filament 
of  strength  n  along  c.  The  strength  of  the  vortex  filament  is  just  the 
local  (edge)  value  of  the  doublet  strength,  i.e., 

u  =  y  (on  c)  (A-2) 

The  vorticity  in  the  sheet  is  a  vector  everywhere  tangent  to  the  curves  of 
constant  y  and  has  an  intensity  equal  to  the  magnitude  of  gradu.  Spec i f i - 
cally,  if  (o  is  the  vector  vortex  strength  on  S,  then 

<o  =  -  n"  x  grad^  y  (A-3) 

Since  y  is  defined  only  on  S,  only  the  tangential  component  of  its 
gradient  is  defined.  However,  it  is  clear  from  the  form  of  (A-3)  that 
the  normal  component  of  the  gradient  does  not  affect  the  result. 


DISCUSSION:  The  Biot-Savart  law  gives  the  velocity  at  (x,  y,  z)  due  to 
a  vortex  filament  of  variable  strength  n  lying  along  any  curve  c  as 


Qds 


(A— 4 ) 


where  s  denotes  arc  length  along  c.  Thus,  the  velocity  due  to  the  vortex 
filament  whose  strength  Is  given  by  (A-2)  and  which  lies  along  a  closed  curve 
c  is 


(A-5) 


The  expression  for  the  velocity  due  to  a  vortex  sheet  Is  obtained  from  (A-4) 
by  writing  the  vector  vortex  strength  u  =  ilt,  so  that  (A-4)  becomes 


■^ds  (A-6) 

FT 

Now  simply  redefine  w  as  a  surface  density  instead  of  a  linear  density  and 
change  (A-6)  to  a  surface  integral  over  S.  This  gives  the  velocity  at 
(x,  y,  z)  due  to  a  vortex  distribution  of  strength  on  S  as 
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(A— 7 ) 


where  dS  is  an  elemental  surface  area  on  S.  For  the  particular  vortex 
strength  given  by  (A-3)  this  becomes 


or 


v 

h) 


(n  x  grad  u) 

- : — tr 


R' 


x  If 


dS 


V 

u ) 


•  f)grad?w 


—  (if  •  gradfp)n 

7  *“"■ 


dS 


(A-8) 


(A-9) 


To  obtain  the  velocity  due  to  the  dipole  sheet,  start  with  the  point 
source  potential 

*s  »  £  (A-l^) 

and  generate  the  dipole  potential 

<fD  =  n“  •  grad^  $$  (A-ll) 

n  . 

=  "7' 

where  n  is  the  unit  vector  along  the  axis  of  the  dipole,  and  in  this 
application  the  axis  is  along  the  normal  vector  to  S.  The  velocity  due  to 
the  dipole  is 


vD  (point)  =  -gradx*D 


grad  (n  •  If)  -  fn  .  "R)  grad 
rj  x 


X 


(A-12) 


The  first  term  above  may  be  evaluated  with  the  help  of  a  standard  vector 
differentiation  formula  taking  advantage  of  the  fact  that  n  is  independent 
of  x,  y,  2  and  the  fact  that  curlxf  =  0.  The  result  is 
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Vp  (point)  =  (?  •  gradx)ft-  (rT  •  I?)  grad/ LA 


1 1  3  2-dJlf 
fT  R5 


(A- 13) 


The  simple  form  of  the  second  form  of  ( A-3 )  is  due  to  the  simple  form  of  if. 
Thus,  the  velocity  at  (x,  y,  z)  due  to  a  normal  dipole  distribution  of 
strength  g  on  S  is 


r  +  3  n— J  udS 


(A-14) 


The  proof  of  the  theorem  consists  of  showing  that  v^  from  (A-14)  equals  the 
sum  of  vr  from  (A-5)  and  -v  from  (A-9).  This  is  done  by  starting  with 
(A-5)  and:  (a)  writing  out  the  line  integral  explicitly  in  terms  of  components, 
(b)  applying  Stoke1 s  theorem  to  each  component  separately  to  obtain  surface 
integrals  over  S,  and  (c)  manipulating  the  result  to  obtain  the  desired 
equality.  The  details  are  somewhat  lengthy.  A  more  concise  proof  should  be 
possible. 


DETAILS  OF  THE  PROOF:  For  a  point  on  the  curve  c  n,  4  are  func¬ 
tions  of  the  arc  length  s  along  the  curve.  The  unit  tangent  vector  to  c 
is 


(A- 15) 


Taking  the  cross  product  with  R  from  (A-l)  and  putting  the  result  in  (A-5) 
gives 

vr  =  1  ^  |  0  dj;  +  (z  —  c)dn  —  (y  —  n)  dc 

T  (j6  |~  U  -  5>dc  +  0  dn+^y(x-Odc  (A-16) 


(y  ~  n)d£,  —  (x  —  c)dn  +0  dc 
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Differentiation  gives 


s 

H 


iT3iL 


+  3 


( A- 1 7 ) 


and  siniiiar  formulas  for  the  n  and  z,  derivatives.  These  are  used  below. 


Stokes  theorem  in  component  form  is 


(A- 18) 


where 


n  =  n11  +  n2wl  + 


(A-19) 


This  theorem  must  be  applied  to  the  components  of  (A-16)  separately.  The 
result  is 


(A-20) 


[+(y-n)(?^+3  Vu 

4 


<k~:i  0) 


)  ni 


(y~")(F^+3jLir")~^Mds 


Certain  terms  can  be  collected  at  once.  The  y/R  terms  add  to  give 


J fa* 


(A-Zl) 


The  coefficient  of  n-ji  Includes  the  term 

2 

-%  C(y  -  n)Z  +  (z  -  c)2]  -  -S  ^  ♦  3ij  (*-22) 

R5  R“  R* 

Similar  terms  occur  In  the  coefficients  of  n2J  and  Separating  those 

terms  and  collecting  gives 


_ 

vr  *  2 


2J^"ds~3-fi?"dS 

+  3  “  5^nl  +  C*  “  0(y  ~n)«2  Mx  -  Ui*  -  c)n3 

+  t  (y  -  n5(x  -  ^)n1  ♦  Cy  -  r;)2n2  +  (y  -  n)(z  -  cin3J 
+  t  (z  -  c)(x  «  +  (z  -  c)(y  -  n)n2  ♦  (z  -  c)2n3||dS 

-  j^4{f[i;x "  ■’  « *i) +  ty  ~ n>  fir  ni +  lz  -  <>  &  ni 

-fx_E)  « "i  !-<>-")  Is  n2- "3  <A-23> 
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>  >  *{y 


&  *».*  (z  -0|fn2 


-  <«  ",  -  j(y  -  ")  ft  »2j  -  (z  -  O  %  «j] 

.f[.N_0  £.3  (<*-«)  {J 

-(x  -  f?  fr  ",  -  (y  -  f*  "2  ~  S'1  -oft  jj. 


SjL  n0 
5;;  2 


(A— 23) 

/  • 

% 

In  the  fourth  Oast)  integral  the  terns  in  dotted  brackets  •  *  have  been  added 
and  subtracted.  I«*  t’?€  third  integral,  if  (x  —  r,)  is  factored  from  the 
first  line,  (j  -  >  )  fro r  the  second  line,  and  {?  -  O  from  the  third  line, 
the  regaining  terms  are  'dentical  in  sli  three  cases v  namely  (rT  •  if).  In 
the  fourth  integral,  the  odd  mesbered  lines  are  identical  except  for  the 
component  of  o,  and  these  three  lines  add  together  to  give  (R  •  grad^u)n. 

In  the  fourth  integral  the  even  ramborsa  11r*es  are  identical  »»xcept  for  the 
differentiation  variable,  and  these  throe  lines  edd  together  to  give 
•(rf  •  *~)  gratify.  Using  all  these  results  (»\-23)  b*cv*es 


gradin’  -  C<f  •  £)grad  u]  dS 


Thus  using  'A-?)  and  (A- 14) 


(A-24) 


(A— 25) 


as  required 
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APPENDIX  B 

LITERATURE  REVIEW  OF  SHAPES  OF  TRAILING  VORTEX  WAKES 

As  part  of  the  present  work  a  literature  search  has  been  conducted  into 
the  problem  of  locating  the  trailing  vortex  sheet.  The  idea  is  that  the  more 
Information  that  can  be  collected  on  this  matter  the  more  accurate  will  be  the 
specification  of  the  wake  to  the  program  and  thus  the  more  accurate  will  be 
the  calculated  pressures.  In  view  of  the  results  of  section  8.5,  it  appears 
that  the  location  of  the  wake  c "  a  lifting  body  is  not  very  important  as  far 
as  the  surface  pressures  on  that  body  are  concerned.  Wake  position  may  be  of 
greater  interest  In  the  case  where  one  body  is  generally  downstream  of  another 
lifting  body. 

It  is  fortunate  that  the  position  of  the  wake  does  not  appear  to  be 
critical,  because  the  literature  has  proved  very  disappointing  in  this  regard. 
First,  there  are  very  few  articles  on  this  subject.  Second,  most  of  those 
few  deal  with  the  asymptotic  wake  location  many  chord  lengths  behind  the  wing. 
This  is  the  Important  region  for  determining  the  effects  of  a  wake  on  another 
aircraft,  but  the  wake  position  at  such  remote  locations  seems  unlikely  to 
affect  the  surface  pressures.  Third,  the  handful  of  articles  that  discuss 
the  wake  in  the  first  few  chord  lengths  behind  the  wing  are  to  some  extent 
contradictory.  Some  of  the  applicable  articles  are  discussed  below. 

Reference  17,  an  experimental  study  of  straight  wings  of  fairly  high 
aspect  ratio  on  a  fuselage,  reports  that  the  wake  vorticity  is  essentially 
all  concentrated  into  the  tip  vortices  right  from  the  beginning.  The  tip 
vortices  separate  from  the  wing  tip  at  about  the  quarter-chord  (not  the 
trailing  edge)  and  go  straight  downstream  parallel  to  the  freestream  direction, 
i.e.,  they  do  not  follow  what  are  normally  thought  of  as  the  streamlines  cf 
the  flow. 

Reference  18  proved  very  encouraging.  The  configuration  was  a  swept  wing 
on  a  fuselage,  and  the  study  was  both  theoretical  and  experimental.  The  wake 
behind  the  wing  was  examined  from  the  trailing  edge  downstream  to  a  distance 
equal  to  one  span.  Various  theoretical  models  were  considered.  One  model 
consisted  of  exactly  the  model  used  in  many  of  the  cases  of  this  report. 
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Specifically,  the  wake  was  taken  to  lie  straight  back  in  the  wing  midplane 
and  the  spanwlse  vortlcity  distribution  was  the  same  as  at  the  wing  trailing 
edge.  Downwash  computed  by  this  model  gave  excellent  agreement  with  experiment 
much  better  than  a  moael  that  considered  the  wake  to  be  rolled  up  into  tip 
vortices. 

Reference  19  presented  the  results  of  numerical  computations  for  wake 
locations  behind  Isolated  wings,  both  straight  and  swept.  The  rolled  up 
portions  of  the  wake  near  the  tips  lay  essentially  straight  back  in  the 
freestream  direction.  The  wake  center  line  lay  much  lower,  but  the  vorticity 
was  quite  weak  in  the  whole  region  near  the  centerline. 

However,  reference  20  contradicted  this  last  result.  Based  on  experi¬ 
mental  studies  of  swept  wings,  the  authors  showed  that  the  wake  centerline 
lay  essentially  straight  back.  This  report  contains  a  large  amount  of  down- 
wash  data  that  is  difficult  to  apply  to  wake-shape  estimation.  This  occurs 
in  many  reports. 
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Figure  26.  Spanwise  distributions  of  section  lift  coefficient  calculated  for  a  swept  tapered 
wing  at  8  degrees  angle  of  attack  using  various  numbers  of  lifting  strips. 


Figure  27.  Spanwise  distributions  of  bound  vorticity  on  a  swept  tapered 
wing  at  8  degrees  angle  of  attack  computed  by  the  two  bound 
vorticity  options. 


Figure  28.  Spanwise  distributions  of  section  lift  coefficient  on  a  swept 
tapered  wing  at  8  degrees  angle  of  attack  computed  by  the  two 
bound  vorticity  options. 
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Figure  30. 
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Spanwise  distributions  of  bouna  vorticity  on  a  swept  tapered 
wing  at  8  degrees  angle  of  attack  computed  with  two  orders 
for  the  input  points. 
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Spanwise  distributions  of  section  lift  coefficient  on  a  swept 
tapered  wing  at  8  degrees  angle  of  attack  computed  with  two 
orders  for  the  input  points. 
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stributions  used  for  the  wall. 


igure  32.  Calculated  effects  of  a  finite  and  an  infinite  plane  wall  on  the  spanwise  distribution  of 
section  lift  coefficient  for  a  wing  of  rectangular  planform  at  10  degrees  angle  of  attack. 


EDGE 


EDGE 


Figure  33.  Two  element  distributions  on  a  wing  of  rectangular  planform 
mounted  on  a  round  fuselage. 


140 


(a) 
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Figure  34.  Comparison  of  results  calculated  for  a  rectangular  wing  mounted 
on  a  round  fuselage  using  two  different  element  distributions 
at  6  degrees  angle  of  attack,  (a)  Spanwise  distributions  of 
section  lift  coefficient,  (b)  Chordwise  pressure  distributions 
at  61.67  percent  semispan. 
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Figure  35.  Geometry  of  an  extreme  test  case  with  a  large  flap  deflection. 
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Figure  38.  Comparison  of  calculated  and  experimental  results  on  a  rectangular  wing  mounted  as  a  midwing 
on  a  round  fuselage  at  6  degrees  angle  of  attack,  (a)  Spanwise  distributions  of  section 
lift  coefficient.  Chordwlse  pressure  distributions  at:  (b)  19.3  percent  semispan, 

(c)  60.0  percent  semispan,  (d)  92.4  percent  semispan. 
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Figure  40.  Comparison  of  calculated  and  experimental  results  on  a  supercritical  wing  mounted  as  a  high 
wing  on  a  fuselage  at  7  degrees  angle  of  attack,  (a)  Spanwise  distributions  of  section  lift 
coefficient.  Chordwise  pressure  distributions  at:  (b)  15  percent  semispan,  (c)  25  percent 
semi  span,  (d)  60  percent  semispan. 


A  conventional  wing  mounted  as  a  low  wing  on  a  fuselage 
(a)  The  complete  configuration,  (b)  Airfoil  section  of 
the  wing. 
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Figure  42.  Comparison  of  calculated  and  experimental  results  on  a  conventional  wing  mounted  as  a  low 
wing  on  a  fuselage  at  6.9  degrees  angle  of  attack  and  a  freestream  Mach  number  of  0.5. 

(a)  Spanwise  distributions  of  section  lift  coefficient.  Chordwise  pressure  distributions 
at:  (b)  15  percent  semispan,  (c)  25  percent  semispan. 


Figure  46.  Comparison  of  calculated  results  on  a  clean  wing,  a  wing  with  tip  tank,  and  a  wing  with  pylon- 
mounted  external  store,  all  in  the  presence  of  a  round  fuselage  at  6  degrees  ang.e  of  attack. 

(a)  Spanwise  distributions  of  section  lift  coefficient.  Chordwise  pressure  distributions  at: 

(b)  92.4  percent  semispan,  (c)  53.3  percent  semispan,  (d)  58.3  percent  semispan. 


Figure  47.  A  rectangular  wing  cf  aspect  ratio  1.4  with  endplates. 
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Figure  48.  Comparison  of  calculated  results  on  a  rectangular  wing  at  10  degrees  angle  of  attack 
with  and  without  end  plates,  (a)  Spanwise  distributions  of  section  lift  coefficient, 
(b)  Chordwise  pressure  distributions  at  the  center  section. 


Figure  49.  A  rectangular  wing  of  aspect  ratio  1.4  in  a  rectangular  wind  tunnel. 
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Figure  50.  Comparison  of  calculated  results  or.  a  rectangular  winq  .  t  ’0  .'.-'.'■r 
with  and  without  the  wind  tunnel  side  walls.  ‘ a )  Spanns*  n1 
lift  coefficient,  (b)  Chordwise  pressure  distributions  ar  the 


Figure  51.  A  rectangular  wing  of  aspect  ratio  1.4  with  endplates  in  a  rectangular 
wind  tunnel. 
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Figure  52.  Comparison  of  calculated  results  for  a  rectangular  wing  at  10  degrees  angle 
of  attack  in  free  air  with  those  for  the  same  wing  with  endplates  in  a  wind 
tunnel,  (a)  Spanwise  distributions  of  section  lift  coefficient,  (b)  Chordwise 
pressure  distribution',  of  the  center  section. 


